{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**This notebook visualizes the relationship between the initial 2D reprojection errors and the final pose errors (Figure 7).**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "from pathlib import Path\n",
    "import pickle\n",
    "import numpy as np\n",
    "from tqdm import tqdm\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "\n",
    "from pixloc.settings import DATA_PATH, LOC_PATH, EVAL_PATH\n",
    "from pixloc.run_Aachen import default_paths\n",
    "from pixloc.localization import Model3D\n",
    "from pixloc.pixlib.geometry import Camera, Pose\n",
    "from pixloc.utils.parser import parse_image_lists\n",
    "from pixloc.utils.eval import cumulative_recall"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "[09/15/2021 12:24:50 pixloc.localization.model3d INFO] Reading COLMAP model /home/psarlin/work/geometry/FeatureMetric-Pose-Estimation/outputs/hloc/Aachen/sfm_superpoint+superglue.\n",
      "[09/15/2021 12:25:02 pixloc.utils.parser INFO] Imported 922 images from Aachen_hloc_superpoint+superglue_netvlad50.txt\n",
      "[09/15/2021 12:25:08 pixloc.utils.parser INFO] Imported 922 images from refinement_Aachen_retrieval-hloc_unet-md-interp75-unc-damp_MS4.1-n3-it100.txt\n",
      "[09/15/2021 12:25:08 pixloc.utils.parser INFO] Imported 98 images from night_time_queries_with_intrinsics.txt\n",
      "[09/15/2021 12:25:08 pixloc.utils.parser INFO] Imported 824 images from day_time_queries_with_intrinsics.txt\n"
     ]
    }
   ],
   "source": [
    "paths = default_paths.add_prefixes(DATA_PATH/'Aachen', LOC_PATH/'Aachen')\n",
    "sfm = Model3D(paths.reference_sfm)\n",
    "\n",
    "hloc_path = paths.dumps / 'Aachen_hloc_superpoint+superglue_netvlad50.txt'\n",
    "hloc_poses = dict(parse_image_lists(hloc_path, with_poses=True))\n",
    "with open(str(hloc_path)+'_logs.pkl', 'rb') as file:\n",
    "    loc_logs = pickle.load(file)\n",
    "    \n",
    "results_path = EVAL_PATH / 'pixloc_Aachen.txt'\n",
    "pred_poses = dict(parse_image_lists(results_path, with_poses=True))\n",
    "with open(str(results_path)+'_logs.pkl', 'rb') as f:\n",
    "    pred_logs = pickle.load(f)\n",
    "    \n",
    "qnames = list(loc_logs['loc'].keys())\n",
    "qdata = parse_image_lists(paths.query_list, with_intrinsics=True)\n",
    "qdata = {Path(q).name: cam for q, cam in qdata}\n",
    "assert set(Path(q).name for q in qnames) == set(qdata)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Get the inlier ranking"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 922/922 [00:09<00:00, 93.06it/s]\n"
     ]
    }
   ],
   "source": [
    "init_poses = {}\n",
    "p3d_ids = {}\n",
    "conf = pred_logs['configuration']['refinement']\n",
    "for qname in tqdm(qnames):\n",
    "    loc = loc_logs['loc'][qname]\n",
    "    dbids = loc['db']\n",
    "    inliers = loc['PnP_ret']['inliers']\n",
    "    \n",
    "    ninl_dbs = sfm.get_db_inliers(loc, dbids, inliers)\n",
    "    dbids = sfm.rerank_and_filter_db_images(dbids, ninl_dbs, conf['num_dbs'], 10)\n",
    "    key = Path(qname).name\n",
    "    if len(dbids) > 0:\n",
    "        init_poses[key] = sfm.dbs[dbids[0]]\n",
    "        p3did_to_dbids = sfm.get_p3did_to_dbids(dbids, loc, inliers)\n",
    "        p3d_ids[key] = list(p3did_to_dbids.keys())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compute initial and final 2D reprojection errors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 919/919 [00:06<00:00, 137.14it/s]\n"
     ]
    }
   ],
   "source": [
    "def project(p3d, im, cam):\n",
    "    T_w2cam = Pose.from_colmap(im)\n",
    "    p2d, mask = cam.world2image(T_w2cam * p3d)\n",
    "    size = pred_logs['configuration']['network']['resize']\n",
    "    max_ = max(cam.size.numpy())\n",
    "    if max_ > size:\n",
    "        p2d = (p2d + 0.5) / max_ * size - 0.5\n",
    "    return p2d.numpy(), mask.numpy()\n",
    "\n",
    "def pose_error(names, poses_gt, poses):\n",
    "    T_w2gt = Pose.stack([Pose.from_colmap(poses_gt[q]) for q in names])\n",
    "    T_w2cam = Pose.stack([Pose.from_colmap(poses[q]) for q in names])\n",
    "    dt = np.linalg.norm((T_w2gt @ T_w2cam.inv()).t, axis=-1)\n",
    "    return dt\n",
    "\n",
    "err2d_init = []\n",
    "err2d_pred = []\n",
    "names = []\n",
    "for q in tqdm(init_poses):\n",
    "    p3d = np.stack([sfm.points3D[i].xyz for i in p3d_ids[q]])\n",
    "    cam = Camera.from_colmap(qdata[q])\n",
    "    p2d_gt, mask = project(p3d, hloc_poses[q], cam)\n",
    "    if np.count_nonzero(mask) == 0:\n",
    "        continue\n",
    "    \n",
    "    p2d_pred, _ = project(p3d, pred_poses[q], cam)\n",
    "    e2d_pred = np.linalg.norm(p2d_gt - p2d_pred, axis=1)[mask].mean()\n",
    "    p2d_init, _ = project(p3d, init_poses[q], cam)\n",
    "    e2d_init = np.linalg.norm(p2d_gt - p2d_init, axis=1)[mask].mean()\n",
    "    \n",
    "    err2d_init.append(e2d_init)\n",
    "    err2d_pred.append(e2d_pred)\n",
    "    names.append(q)\n",
    "\n",
    "err2d_init = np.array(err2d_init)\n",
    "err2d_pred = np.array(err2d_pred)\n",
    "err_init = pose_error(names, hloc_poses, init_poses)\n",
    "err_pred = pose_error(names, hloc_poses, pred_poses)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Define bins and count"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1.0,\n",
       " 0.8461538461538461,\n",
       " 0.9574468085106383,\n",
       " 0.9333333333333333,\n",
       " 0.8933333333333333,\n",
       " 0.8176795580110497,\n",
       " 0.8277777777777777,\n",
       " 0.6752136752136753,\n",
       " 0.410958904109589,\n",
       " 0.35,\n",
       " 0.0,\n",
       " 0.0]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bins = np.logspace(np.log10(20),np.log10(1500), 13)\n",
    "prob = []\n",
    "for i in range(len(bins)-1):\n",
    "    select = err_pred[(err2d_init>bins[i]) & ((err2d_init<=bins[i+1]))]\n",
    "    p = np.mean(select <= 1)\n",
    "    prob.append(p)\n",
    "prob"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Histogram"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Recall of oracle @1m: 4.13%\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdIAAADdCAYAAAAPQfNXAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABEk0lEQVR4nO2dZ5hT1daA3zVDb4KFIgLRiwVsIKgoCuhF1C/2rtiwYFf0okZBBa5IVK5iF7gXQQR7JyiKir0hNsSGGFBAUECl1/X92CeQyWQmJ5NkMmW9z3OeOWfXdU7OZGXvvfZaoqoYhmEYhlE2CvItgGEYhmFUZkyRGoZhGEYGmCI1DMMwjAwwRWoYhmEYGWCK1DAMwzAywBSpYRiGYWRAjXwLUF4UFBRo3bp1y1RXVRGRLEtUvn1m2l5Z6qdTx0/ZVGUyza8MlPc92HuY/TKV7T1ctWqVqqoNukpDVavFUa9ePS0rP//8c5nrVpQ+M22vLPXTqeOnbKoymeZXBsr7Huw9zH6ZyvYeAiu1tO9X6KbwksJ8BVU4NyFfFAYpLFBYrTBNYfeEMk0Uxiv85R3jFRrH5QcU3lFY6f0NJNR/VqFvqXLm8PD1KyMQinQPhCL7x12fGwhF3guEIiMDoUiDnGn5ysqECRAIQEGB+zthQr4lMgzDyBUNgJnAVcDqJPnXAf8CrgD2BRYDryPSMK7MRGAf4Ajv2AcYH5f/H2A+0AFYCAzfnCNyHLAdMDoL91Im/A7XRwDNAQKhyK7ASOAr4ADgzpxIVlmZMAH69oW5c91vs7lz3bUpU8MwqiKqk1G9EdVngE1F8twcdj8gjOqzqM4EzgEaAmd4ZdrhlGdfVD9E9UPgIuAoRHb1WmoHjEP1R2Csdw0ijXBKti+aPzd9fhVpW+Br7/xE4PVoOHgpcCFwdC4Eq7QMGACrVhVNW7XKpRuGYVQvdsQNwl7bnKK6GngHONBLOQBYAXwQV+99YGVcmS+BnogUAL1wAzmAMDAW1e9yJL8v/BobbQIKvfN/As97578B22RbqFyw9dZbM23atDLVXbduHdFo1FfZ7vPmkcyMQOfN4+00+k+nz/Joryz106njp2yqMpnmVwbK+x7sPcx+mUr4HtYQkelx16NUdZTPus29v4sS0hcBLePK/F5kRKmqiCyOq98fNxMaxSnRixA5EDgYOASRCTiF/BFwMap/+5QvK/hVpJ8CNwVCkddxgvf10gO4+eoKz9KlS+nRo0eZ6kajUQKBgL/CrVu76dwEpHnztPpPq89yaK8s9dOp46dsqjKZ5lcGyvse7D3MfplK+B5uUNXOeZVAdT5w1OZrkVrAq7gp4BCwAdgFN+17E3BteYrnd2q3H26R935gaDQc/MlLPxn4MPtiVWKGDoV69Yqn//03fPpp+ctjGIaRP37z/jZLSG8Wl/cbsF2RPUHuvGlcmURCwDuofgAcCjyF6gbgce+6XEk5Ig2EIgW4qd0Do+HgioTs/sDGXAhWaend2/0dMADmzXMj1Kuugvvug0MPhZdfhjKOjPPChAkwYABtYvcydOiWe6wGTJ2VOCOVfXq2T/yOMYwqw884ZXgYbmYTROrgZjZjo8YPcZa/B7BlnfQAoD5F103x6u+GM1jq6KUUADW981psWYYsN/xM7SrwBdAemB2fEQ0H1+RApspP797Flc0pp0CvXnDEEfD003B0JbDRilkgr1rl1n1jFshQrZRprklHWS9evJzZq8qm3E1hGzlBpAHOIBWcUmuNSAdgKarzEBkB3IjId8APwECccdFEAFS/ReRVYCQisWXDkcAkVL9P6EuAUcA1ceug7wEXI/ItcIl3Xa6knNqNhoMKfI/bp2OUlZYt4Z13YK+94Pjj4bHH8i1RakqyQL7xxvzIYxhGRaQz8Ll31AUGe+dDvPw7gLuBB4DpQAugF6rL49o4A2eZO8U7vgTOStJXX5xh0otxaYNwA77puNnTQVm4p7Twa2x0HTA8EIpcBnzpKVcjXbbZBt54A449Fs46C/78Ey6/PN9Slcy8eSWn77cf7Lwz7LKL+xs7GjcuVxENw8gzqtMg6WaFWL7ilNugUsosA8700ddI3Gg1Pu0P4MiUdXOIX0X6FFAH+AzYEAhF1sZnRsPBRtkWrMrSsCFMngynnQZXXOGU6YABUJF8b6rCyJEl5zdoAFttBe+/D48/7srH2G67zcp1q+22cwp3l12gbdvkRlgAEyaww/XXw4IF1XId1jCMyo1fRVqBh02VkDp14Jln4Lzz4KabYNkyGD68YijTX36B88+H11+H3XeHOXNgdZzXr3r14OGHtyi6NWvgp5/gxx/hhx+2/J0yhSYLE3ZG7bBD8VHs99/DoEHUiE0h2zqsYRiVDF+KNBoOjsu1INWOGjVg7Fg3FXrXXW5kOnKkS88HqvDII3D11bBxIzz0EFx0EUycCAMGoPPmIclGi3XqOIW7++7Fmpw7cyZt1q8vqmB//NH9iFiypGRZYuuwpkgNw6gE+P7WDoQitYHeOOtdBb4BHo+Gg2tLrWiUTEEB3HMPNGkCQ4Y4ZTpxItSuXb5yLFjgRoGRCHTvDmPGwE47uTzPAnluGTaRa4MGzml/x47FM5cudUq1S5fklefNc1bO3bq5Y7/90urbMAyjvPAb/aU98CNwF7A/0AXnyP6HQCjSLmfSVQdEYPBguPtueO45ty1m5cry6VvVbXHZYw94800YMcL9jSnRXLL11rD//tCmTfL8hg1h0SI39d29O2y1Fc1POQUGDoTXXoMViVuaDcMw8oNfz0b34MyZW0fDwYOj4eDBQGucifKIHMlWvejXz02tvvEGHHYYBX/9ldv+Fi2CE0+EM8+E3XaDL75wjiMKyjl+bzJPUPXquanlL790U8AvvghXXIGsXQvhMBx+uJsS339/uPZa5+Ri2bLyldswDMPD77dmV+DGaDi42RGwdz4AOCgXglVLzj3XOWv47DOan3oq/FaSd6wMefppNwqNROCOO+Ddd50BUD7o3RtGjWJDy5ZudN6mDYwatWV9dOut4ZhjYPhwFr74olOYU6ZAKOSmwO+91+Vvsw3bH3mks4R++mn3QyGB+i+8YHFiDcPIOn7XSNcAjZOkb+XlGdnihBMgEqHGscfCQQfB1KnuSz8bLFkCl10GTz4JnTvDuHHQvn122s6E3r35tWtXf2uwDRu6tdNevdz16tXwySfwzjtsnDLFre/ef7/L23XXLWusS5eyzY03brFANutgwzCyhF9F+jIwOhCKXIgLUwPOF+JI4KVcCFat6dmTRY89Rovzz4euXd1WlAwVXt3XX3frjUuXwr//7UZ0+bIQziZ167o11O7dWXTWWQRatoQZM5wXqXfegaeegtGjgSTTL7E4saZIDcPIAL9Tu1fhjI3exY1A1wBv4/wm9suJZNWctR07wttvw6ZNbkRV1sgxy5bB2WfTrG9faN7ctTNwYNVQosmoWbPo2umSJfD55yWXL8l7k2EYhk9SKlIv+sv2uK0vuwAneMeu0XDw+Gg4mGOrmGrMnnvCe+9Bo0Yucky6gclffdWthU6cyJ9XXOGmQPfeOyeiVlgKC6FDh5KtgwsL3Yh1ja1QGIZRNvyMSGPRX5pHw8HZ0XDwZe+YnaKekQ3+8Q9nDNS6tYsc8/LLqev8/TdceCEceaSzbv3oI/685hqoVSvn4lZYhg5lU926RdNq1XLBBPr2devQt91m1r+GYaSNRX+pDKQTOeaNN9xIdswYuO46+OwzZ1hU3endmyW33eZGpjHr4DFj4Oef3TPr0MGtl7Zq5bw72ZSvYRg+8btGGov+0iEQilQAh7DVkFjkmG7dXOSYc84pupVjzBhnkduzp9sW8t57cPvtzoWfAcDK446DaNStO0ejzshIxE2bv/qq20t7/PEuCPtOO8GZZ9Lgu2/yK7RhGBUei/5SmYhFjunaFR59dEv63LlwwQXOU1G/fsmdHBip2XtvGD/ePb8RI2D0aLpMmMCSA3sw97xLWdrl4IoRWMAwjOwgUhtnA1QXF+f097I0Y9FfKht16sAffxRPV4VmzZyrQSMzWrd2gQRuuonZQ4bTavxo9rngFP5utydz+1zK4sOPRquq1bNhVHVEGuJin54O7AfUxMVTVUTm4wKLj0LV91aJlN8GgVCkJvB/OM9GP5VFbiPL/PJL8vTFi8tXjqpOkyZEL7ySeWf3pfnLz9Jm7EPsed0lrL7nNuadfRHzTzidTfXql7tYzSY9S9sRw6jz23zWNG/J7H43sOioE8tdDsOodIhcg/PINwfnA2EosABYDWwN7AEcDLyOyEfAFaj+mKpZP8ZG64FeOOtdoyLQunV66UZGbKpdhwUn9ebDl97hy/vGsrZpC3YdNpCD/9mJne69nZpLyjQbVCbaTI3Q/pb+1F34K6JK3YW/0v6W/jSb9Gy5yWAYlZguQHdU90X136hOQfVrVGej+gmqY1DtAzTDKdrufhr1a2z0HG7vqFERKMnR+9Ch+ZGnulBQwO+HHsH0x17i08deZtm+B7DjqBEc1LMzuw2+jrpz59Bs0rN07dmZf+7Rgq49O5ddwalSuHIltX9bQP0fv2WrGZ+wzduv0+m+2ylcs7pI0cI1q2k7YlgWbtAwqjiqp6A600e5tag+iOp//TTrd6FnHjAwEIocDEwHisT5ioaDd6VqIBCKdAP6A51wi7t9ouHg2Lj8scA5CdU+joaDXeLK1AaG4+a26wJvAJdGw8Fffd5H1SDm0m7AALdNI1nAbSOn/NVxX77q+Aj1fp5N67EPs/3zT9DyqUehoBDZtBHAjRZv/hd15//CXx06U2PFcmr8/Rc1lv9NjRV/U/Pvv6mx/K9i6TX+dn8LNm70LU+d3+bn6lYNw0iBX0V6LrAM2Ms74lFcnNJUNABmAo96RzKmAmfFXa9LyB8BHItTpEu8ficFQpFO0XDQ/7dOVcALuG3kl1U7tuW7wcOZc/l1HHBUV2quWF4kv3DtGtreG05ad0P9Bmxo2IgNDRuxvuFWrG3ajJX/2GVzWix9Q6NGbGjgrve4/BzqJZlKXtN8+5zcn2FUaUSa4aKbNSVxhlb1Qb/N+FKk0XBwx3RkK6GNycBk2Dz6TMbaaDiYNHZYIBTZCjgfN5J93Us7C5gL9MRZWhlGXli3XVNqrEwebFyBGWOeYUOjrVjf0CnFjQ0alsny9/OLr6HLf4YUm97dWKcuhSuWs7FBw7KIbxjVD5Ezgf/iLHaXUdQOSIHsKtJ4AqFIM+D3aDi4Kd26PjgoEIosBv7EOcUfEA0HY6aonXBmyq/FCkfDwV8Coci3wIEkUaQi0hfoC1CrOrvHq6JMnVU05ujixcuZvap4HFK/+ZmypnlL6i4svsqwpsUOLNs/O2F75/YM0qhRoyJWu79368kOzz5Gpz4n8vnDE1i/jTkhMwwfDAXuAIaguiGThnwpUm8LzFDgEtza5C7AnEAocjswNxoO+tbcpfAqzqjpZyAA3Aq86U3brgWaAxuBxE2Ui7y8YqjqKGAUQP369c3q2Mgps/vdQPtb+hcZLW6sU5fZ/W7Iaj+Ljjqx2HaXJd3/yV5XX0jns49jxugnWbv9Dlnt0zCqII2AsZkqUfBvtXsLcDRuE2u8V6NPcOunGRMNB5+IhoMvRcPBr6Ph4MvAkcCuQDAb7RtGrll01InMGjyc1S12QEVY3WIHZg0eXi57PJd0P4wZo5+k1pLf2fesY6j30w8579MwKjkTyJJ+8Tu1ezpwXjQcfDsQisRP6c7EjU6zTjQcXBAIRX4FdvaSfgMKgW2BeGuLZrg4qYaRd5KNFsuLvzrtz2djn6dj39PofPZxfD5yIsv36JAXWQyjEnAN8AIi/wS+BtYXyVUd4rchvyPS7XFGPYnUoAzrrH4IhCLbAi2BhV7SZ7gbPSyuzA5AO+CDXMhgGJWNFbvtzvTxL7Gxfn069TmRJh+9l2+RDKOichFwBM7G5njg5LjjpHQa8qsEvwG6AdGE9FNwCi4lgVCkAdDWuywAWgdCkQ7AUu8YBDyLU5wBYBiwGHgeIBoO/hUIRf4H3OEZJMW2v3yF2zZjGAawus2OTB//Eh37nkbHi8/g6+EP83vP/8u3WIZR0bgJ+BeqGTso9zsiHQzcFwhFBuCmV08OhCKPACHg3z7b6Ax87h11vTY/B4bgjIj2BF4EfgDG4WKgHhANB+M35vXDKdYngfeBFcDR1W4PqWGkYG2zFkwf9wLL2+3BXldfwPbPTcy3SEZVRSSKiCY5Il7+oCR5vyW00R+Rxd7xr4S8joh8h0jdLEteiHMDmDF+95G+HAhFTgFuBDbhjI9m4JSYr9FgNBychtuvUxKH+2hjLXCFdxiGUQobGjdhxn+fZq9+59H+pmug7ib4179SVzSM9NgXp5RitMDNVD4Vl/Y90CPuesvgR2Qv3IDqKJyOmITIa6h+jUghMBq4HNWim6cz5xGgt9d3Rvhe34yGg1MwpweGUanYWL8+XzzwKHuELqdZ//4uBN9tt1lcVSN7JMbwFDkf+JuiinQDqkmd7QC7AV+h+qZX/ysv7WvcLORMVHOxfFcPuACRw3FLhInGRlf6bUhUq8f2ylatWun48ePLVHfdunXl7tAh231m2l5Z6qdTx0/ZxDLL1xTd/rV+w3pq1qhZYv1U+ZWBMt/Dxo10evh+tp80iQVHHcUP/fpBYWHKavYeZr9MPr5PMuGQQw5Zh1NqMUZ5e/SLIyLAT8BkVC/30gYB1+Ec7awFPgZuRHWOl98Ot1TXATci/QJnALQGeBPojOqSbN6T1+9bpearHuK7qeqiSOvXr68rV65MXTAJ0WiUQCCQXYHKuc9M2ytL/XTq+CmbWKa4Z6PFNG3atMT6qfIrA5ncQ892TV2gg2HD4OSTYfx4qF271Dr2Hma/TD6+TzJBRFapqr/AuyK9cDOXHVD90ks7EmgIfIfzaTsQN+LcfbOCFLkYuNpr5W5UH0bkVWA8LkjKEJxNz0BUX8jKjWWRnGxdMQyjAiLipnW33hquvRb+/BOeew4aNMi3ZEbV4ULg081KFED1lSIlXMDsObhoX3d5ZR4GHo4rE4vI8TrOAPUAnCJ9H5FdUF1MrhBpBQxG9Ty/Vfxa7RqGUVXo3x/GjIE33oDDDoOlS/MtkVEVEGmKi841utRyqitwWyp3Tpovsg3ORezFuEDcP6L6LarfAD8C+2dP6KRsTfGQnqViI1LDqI706QONG8Npp0G3bjBlCrRsmW+pjMrNubg10MdLLSVSBze1W9Ia5V3AfahGEemAC1YSoxZFLYTTR+TsFCVap9tkiYo0EIpc47cRP4G9DcOoYBx/PLzyChx7LBx0ELz+OrRtm7qeYSTijIwuAJ7wRpzxecOBl4F5uDXSm4D6OH8Bie30BNoDsWnVT4FdETkGZ4i0K87HeyaMBVZRNGxaPGnP1JY2IvW7V9NvYG/DMCoahx4Kb74JRx7plOmUKbD33vmWyqh89MBN1Z6ZJG8H3Cg15if9I6ALqkXdzjqHCw8Ap6Hq9pmqzvcMkR7CKdKLUF2QoawLgCtRfS5prhsF+/LYF6NERZqNYN6GYVQC9t0X3n0XevWC7t1h0iSnVA3DL6pvUZLDHdXTfLaxGjfiTEwfR7LRa9n5DNgHF7YzqSSU7jyoGGZsZBgGtGsH778PzZo5hTp5cr4lMoxcMRy3b7UkZgO+95BCGsZGgVCkCS5GaGvcgu9mouFgxi6WDMPIM61bu5HpkUe6ddMLL2SHl16CBQtc3tCh0Lt36nYMoyKjWnrYTdWVwNvpNOlLkQZCkS5ABGeRtR0wH+dPcS0uIowpUsOoCjRtCm+9BfvvDw89tOULYu5c6NvXnZsyNYwi+J3avRMXTbwlzm3TobiR6XTg9tyIZhhGXmjUCJJ5AVu1ynlGMgyjCH4V6V7A/dFwUHFe+2tHw8FFwPW4OKKGYVQlfv01efq8eeUrh2FUAvwq0nVx54uANt75CmD7rEpkGEb+aV3CnvSS0g2jGuNXkc7AxZwDmAbcGghFzgHuxYWfMQyjKjF0KNSrVzStdm2XbhhGEfwq0gG4TazgPPf/DtwHNAH65kAuwzDySe/eMGoUG1q2dM7uCwuds/uTTsq3ZIZR4fBltRsNB6fHnf+O2wZjGEZVpndvfu3a1YX8mjwZgkG4804YODDfkhlGhcIcMhiGkZr/+z8Xw/TWW+HHH/MtjWFUKEpzWv8V0D0aDi4LhCJfU7KDX6Lh4F65EM4wjArEiBHOF+8llzgH95KWFzXDqDiIbAXsldQ5g0hXYBaqy/w2V9rU7rM4hwux8xIVqWEY1YDtt3eBwS+/HCZONMcMRmVmE/AKIoejusVdoMjewJs4nwm+Kc1p/eC480Fpi2kYRoVi6qxFaddZvHg5s1fF1et2PPvu+T/qXnkVH+zYiQ2NmxQp37N9s0zFNIzco7ockReBsynqd/csYAqqf6TTnK810kAo8mYgFGmcJL1RIBR5M50ODcOoxBQW8u0td1Ljrz/Z+a5b8y2NYWTCo8DJiDjf8SIFwBm4eKVp4dfYqAcJjuo96gAHp9upYRiVlxXt9mDe2X1p+ewEtvrs43yLYxhl5XVgNXCUd/1PnJ57Od2GSt3+EghF9om73CsQiiyNuy4EDsc5sDcMoxox59L+NHv1JdoNvpaPn5mK1kr2O9swKjCqmxB5DDe9+xxuWvdJVNen21SqfaTTcUZGCryWJH81cEW6nRqGUbnZVK8+3w8cRofLzqbN2IeI9r0q3yIZRll4FPgMkdbA8bhRadqkUqQ74iKFzwH2w3k0irEOWBwNBzeWpWPDMCo3f/ToxaLDguz48N0sOuJYVrcO5Fskw0gP1W8QmYmLbvYrqp+UpZlSFWk0HJzrnZrjBsMwivHDDbeyzQdvs9u/Q3w+6vF8i2MYZeFRYATOFW6ZKM0hwwnAy9FwcL13XiLRcPC5sgpgGEblZW2zFvx0VYhdbxtIs8kvwO4X51skw0iXx3B+4x8pawOljUifAZoDi73zklCc4ZFhGNWQX07rQ/MXn2aX22+CC06FJk1SVzKMioLqUmBwynKlUJpDhoJk54ZhGEUoLOS7wcPZ75TD4YYb4OGH8y2RYZQrpiANw8iY5e32ZN6ZF8LIkfDBB/kWxzDKFV9h1AACocgOQDegKQkKOBoO3uWjfjegP9AJ2B7oEw0Hx8blC3ALLr5pE+Bj4LJoOPhNXJkmuGDix3hJLwFXRMPBP/3eh2EYuWHO5dfR5q3JcNFFMGMG1KyZb5EMo1zw6yKwN/ATMBroh9s7Gjsu99lXA2AmcBVu/2ki1wH/8trcF7c2+3ogFGkYV2YisA9whHfsA4z32b9hGDlkY/36cP/9MHMm3JXyt7Vh5BeRmV4UmIzxOyIdAvwHuKms+0aj4eBkYDJAIBQZG5/njUb7AeFoOPisl3YOTpmeAYwMhCLtcMrzoGg4+KFX5iLg3UAosms0HPy+LHIZhpFFjjkGjjsOBg928Ut32infEhlGSbQHahdLdcr1NlQv89uQ3zXSZsB/c+h8YUechfBm70nRcHA18A5woJd0ALACiF+AeR9YGVemCCLSV0Smi8j0DRs25EJuwzASufdeKCyEyy4DteiLRgVD5BVEBuN2nLRKUqIecFE6TfpVpJOB/dNpOE2ae38T4zwtistrDvweDQc3/2d654vjyhRBVUepamdV7Vyjhu/lYMMwMqFVK7j1Vnj1VXj66XxLY+QakUGIaMLxW1y+eGUWILIakWmI7B6XXxuR8Yj8jcgPiPRMaP8KRCZmUeKvge44r32fIPInIm8jMgKR84CrgYXpNOhXu7wO3B4IRXb3hCji1NccMhiGUYTLL4fx4+Gqq6BXL2jcON8SGbnle1yUsBjxs5cx+5dzvXI3A68jsiuqy3EGpp1ws45HAhMRaYaqItLKq7tf1iRVvQ4AkbVen9sDHbwjiNOL16XTpF9FOtL7e2MyscjcIUPs10szYF5cerO4vN+A7QKhiMRGpd7aatO4MoZhVAQKC91WmP32gxtvhAcfzLdERm7ZgGrx72GRzfYvqD7rpRWxfwHaAS95fm/nAHcC2+J8uz8IDEJ1cQ5kbuBFepkBTEqQu2fSGiXgS5GWg0OGn3HK8DDgU4BAKBKLdXqtV+ZDnOXvAWxZJz0AqE/RdVPDMCoCnTrBFVe4NdOzz4YuXfItkZE7dkJkAbAWt3XxRlTnkMT+BdXViMTsX0YCXwJnIVIXF5pzIfAHIqcAdVEdmxOJE8OlibQE+nhHG9LYHipaTsYAgVCkAdDWu/wACOP2gS6NhoPzAqHI9bgRbx/gB2Agbt/qrtFwcLnXxivADripAIBRQDQaDh6dqv9WrVrp+PFl2ymzbt06apVzvMVs95lpe2Wpn04dP2UTyyxfU9SAbP2G9dSsUfLexVT5lYHyvod0+2tYp+h3T+GqVex3zjmsb9SIz0aOZO2mTVXuPUy3TD6+TzLhkEMOWYdb0osxSlVHbb4SORJoCHyHmyEcCOwG7A7sijMKbYPqvLg6Y4CWqB6OSE2c0/j/A/7ArVHOxI0UewG9veM3oC+q32Xt5kQKgWOB83FK/FfgbeBMVH3PtJbmtP4a4MFoOLjGOy8RPw4ZgM7AW3HXg71jHG7u/A6gLvAAWxwy9IopUY8zgPuAKd71S/jcx7p06VJ69Ojhp2gxotEogUCgTHXLSrb7zLS9stRPp46fsollps4qapu2fPFiGjVtWmL9VPmVgfK+h3T769G+WfHE0aOpffzxdP/8c6Inn1zl3sN0y+Tj+yRDNqhq5xJzVV8pci3yES705jnARylbdyPDoltNREbjBkq7Aafi1lBPx/kN2DcN2ZMjsitOeZ4N1AKeBg5B9V1E9gDOTKe50oauV+CU3BpKD96tQEpFGg0Hp+GspErKV2CQd5RUZhlp3qBhGHnmuOPc/tJBg6hxwAFQuZSIkS6qKxD5BtgZeMFLLc3+pSgi3XEDr0twA6wIqssRmQCMRKShZ6RUNkTexSnmCHAxMBnVdfF3kG6TpTmt3zHZuWEYRtrcdx+0b8/WN98Mb7wBUuJvaqOyI1IHN5J8iyT2L15+vP1LfN3awEPAOahuQKSALds0Y/PhmRq3dgWeBO5H9f0M2wLSWEw1DL/EplwXL17O7FWJW4OT46dsOu0ZFYzWreHf/6beNdfAs8/CSSflWyIjW4gMB17GjTibAjfhjEDHeVtYRgA3IvIdW+xfVuBcviZyEzAF1U+96/eAuxEZC5wCfIPqnxlK3Am4AJiEyDJPjgmoflvWBi36i2EY5cMVV7C2fXu48kr46698S2Nkjx2Ax3F7RJ/DWe52QXWul38HcDfO/mU60ALoVWx61q1NnopTtDGe8463cMFKzslYWtXPPfd/LXCBUroCMxGZgcjVuH2laWGK1DCM8qFGDZbcdhv89hsMHJi6vFE5UD0N1e1RrYVqS1RPRHVWXL6iOgjVFqjWQbU7qjOTtDMT1Z1RXRmXtgnVq1BtguoeqH6WRbnXoDoe1UNw1sWv4aabp5ResTimSA3DKDfW7b2383r0wAPwySf5FscwHKqzUQ3hfO8eT6KDhhSUqEgDocjNgVCknnfe2vMiZBiGkRm33gotWri4pRZMwqhIqG5E9UVUj02nWmkj0ptxnoTAWV5tV1bZDMMwNtOokfN29MUX7q9hlBci/negOGf7yaLDFKM0q935wEmBUCSC2/+5g+e2rxjRcHBesnTDMIyknHACHHUUhEIuCPiCBc6yd+hQ6N0739IZVZcPEYkA/0X1w6QlRJoApwFX4gyk7k/VaGmKdKjXwH24DaqfJikjZMdpvWEY1QkR6NkTJk2C+fNd2ty50Nfz/mnK1MgNuwEDgAgim4DPgAU4x0NNcMG+2wGfAP1Q9WV4VOLUbjQcHIXbE9QJpzCPxIWyiT/2JZvhbQzDqD7cfXfxtFWrYMCA8pfFqB6o/onqtUBLnFejb4HGOOf6G3De/Dqi2tWvEoUUDhmi4eCfwBeBUKQP8HY0HFxbNukNwzASmFfCilBJ6YaRLVRXA894R8b4DaM2DiAQihyKG/oqMCsaDr5VakXDMIySaN3aTecmogpdu8I558App1hQcKPC42sfaSAUaRkIRT4BXgeuB0LA1EAo8nEgFEnbC4RhGAZDh0K9ekXT6taFU0+FZcvc9pjmzeG002DyZNsqY1RY/DpkuBfYCLSNhoOtouFgK5xn/41enmEYRnr07g2jRkGbNs74qE0bGD0anngCvvkGPv0ULrwQpk6FYBBataLJ0KHw1Vf5ltwwiuBXkR4GXBYNB3+OJUTDwTk48+DDciGYYRjVgN69IRqFTZvc35i1rgh07uyixixYAM8/D1260GjcONh7b+jY0RkrLbIgBkb+ScdFYLIYbWnHbTMMw0iLWrVcTNPnn+eXjz5yyrVGDbjmGmjZEo4+Gp55BtasybekRjXFbxi1N4D7AqHI6dFw8BdwbgOBEV6eYRjG5hB6JZFpKLzFf9ek6aEn0/Pyy2HWLHj0URg/3u1HbdzYraeefTZ06WIxT43iiLzku6zqMX6L+h2RXomLLzcnEIrMDYQic4GfvLQrfQtmGIaRLdq3h3DYbZd57TW3jjpuHBx4IOy2mzNmmjuX+i+8AIEAFBS4vxMm5FlwI48sSePwjd/tL78EQpF9gJ44zxAA30bDwanpdGYYhpF1CgvhsMPc8fffLnD4uHEuVNvAgWxbUODWYMG8J1V3VPvkolm/U7tEw0HFbX95PReCGIZhZEyjRtCnjzt+/hk6dkQSg4jHvCeZIjWyhG9FahiGUanYcUc3Qk2GeU+qnuRojdQUqWEYVZeSvCe1bl3+shgVgbTWPv1iitQwjKrL0KFsuvBCClav3pJWs6YzRDKqH/leIzUMw6h09O7Nkt9/Z7sRI9x0bu3abl/qCSfkWzKjIiBSAxfBrDVQKy5HUR3vt5m0FWkgFGlMwraZaDi4NN12DMMwyoOVxx3Hdv36uYu33oJDD4WxY+GSS/IplpFvRHYDXsaFUBOcy9sawHpgLZBdRRoIRdoADwM9KKq1LbC3YRiVhx49YP/94c47nR/fGjYpV40ZgQvs3QH4zfu7FfAQMDCdhvy+RY/ggp+ej4smbq4BDcOofIhAKATHHw9PPQVnnJFviYz8sS/QHdWViGwCaqA6A5HrgPuAvfw25FeR7gd0iYaDM9OX1TAMowJxzDFbvCKdfrq5Eqy+CLDKO/8daAl8D/wKtE2nIb8uAn8GaqfTsGEYRoWkoACuvx6+/trFOTWqKzOBvb3zT4DrEekODAZmp9OQX0V6FTAsEIqkpaUNwzAqJKef7vaShsP5lsTIH0Nxo1Jwa6KtgbeAXqTpQ97v1O6LuBHp94FQZC1QJFR9NBxslE6nhmEYeaVmTejfH668Et57Dw46KN8SGeWN6pS48zlAO0S2BpahmpYdkF9Fenk6jRqGYVR4zj8fhgyBYcMgEsm3NEZFQLVMWzn9Rn8ZV5bG0yEQigwCbklIXhQNB5t7+eLl9wWaAB8Dl0XDwW9yLZthGFWQevXgqqvgppvgq69gL99GmkY8IjcAJwC74vZffgTcgOrMuDJjgXMSan6Mape4MncB5wIrgRCqE+LyjgauBw5Od7RYitxDgV9QfTgh/WKgJao3+W3K9yaqQChSG+gNtMdtf/kGeDwaDq7124YPvsftVY2xMe78OuBfuAf9PXAz8HogFNk1Gg4uz6IMhmFUFy67DG6/3a2VTpyYb2kqKz2AB4FPcWuOQ4CpiLRPGOFNBc6Ku163+cwpyjNw65M7A2MQmYLqH4g0BO4GjsmaEnWcBZycJP0z4AYgu4o0EIq0B14FGgFfe8kXAoMDocgR0XDwW78dpmBDNBz8LUn/AvQDwtFw8Fkv7RxgMe7hj8xS/4ZhVCeaNIGLL4a77oJbb3UWvUZ6qB5e5FrkLOAvoCvOc1CMtagW+373aAdMQ3U6MB2RETiPQ38AtwGPoTory5I3xW17SWQJ0CydhvyOSO8BPgfOioaDfwMEQpFGwGM47xCHl1w1LXYKhCILcNMDHwM3RsPBObgH2hx4LVYwGg6uDoQi7wAHUoIiFZG+uKlgWrZsybRp08ok1Lp164hGo2WqW1ay3Wem7aVTf/UaZ4tWZ8N6Vq9Z5KuOn7KpymSaXxko73vIdn+ZtherP22x/9/uqd7dWvvvT5fCQhZefTU/XHZZyvfcz/9CaWXy8X2SITVEZHrc9ShVHVVK+Ya4HSHLEtIPQmQx8CfwNjAA1cVe3pdAX0SaADsBdYHZiHQBDgH2yfw2ijEPOBiYk5DeDbeX1Dd+FWlXYN+YEgWIhoN/B0KRAbj58GzwMW7a9jvcL4WBwAeBUGR3nBIFSPwPXITbRJsU78MeBVC/fn3t0aNHmQSLRqMEAoEy1S0r2e4z0/bSqT91lvuYli9eTKOmTX3V8VM2VZlM8ysD5X0P2e4v0/Zi9Xu09z9g8PXuvvYaLR99lE0DB9Jq330zbq+0Mvn4PsmQDaraOY3y9wBfAB/Gpb0KPIfzSRAAbgXeRKQTqmtRnYLIY7jp4dW49dQVuO/vi4E+iPTDOVC4AtUPMrojx0jgbkRqAW96af8EhgG3p9OQX0W6BuciMJGtvLyMiYaDr8RfB0KRj3C/FM4he8raMAyjONddB//7H43GjIEUitQoBWcwdBBwEKpbbFxUn4gr9TUinwFzgSBOwYLqIGBQXFsDgA9w08RDcL5w9wSeQmQnVNeRCar/QWRb4F62+JBfB9yD6h3pNOVXkb4MjA6EIheyRakdgNPo/iOOp0E0HFwRCEW+wS08v+AlN8MNx4m7LmnO3UhCbLSYLosXL2f2qso9LWoYJdK2LZx0Eg0fe8xth2ncON8SVT5E7gZOAw7x9mWWjOoCRH7Ffb8na2sXnG/3jrjB1DuoLgQWIlIbZyH8ddK66aB6AyK34oxoAb5FdUW6zaTj2ehH4F3cCHQNbo77B5wRUNYJhCJ1gN2AhbjpgN+AwxLyD8b9YjEMw8iMUIiCFSvgoYfyLUnlQ+Qe4HTgUFS/81F+W9yy3MIkeYKb0u2P6l84PVUzLq8m2Yw4proS1U+9I20lCv73kf4JHBsIRXbGKTeAb6PhYFr+CEsjEIoMx4185+HWSG8C6gPjouGgBkKREcCNgVDkO5wCH4ibQzebdcMwMqdjR1Z17069ESOgXz+oWzffElUORB7AbSU5DliGSMymZQWqKxBpgJuyfRanOAO4dcjFwPNJWjwf513oOe/6PWAIIgfhfOOux22BzIbsR+IcDu0IHI7qL4hcAPyM6ht+m0krGF80HPwRNzLNBTsAjwPb4kySP8JFnJnr5d+Bs+R6gC0OGXrZHlLDMLLFXxdfTL3TT4dHHoFLL823OJWF2INKVDyDcQp0I25t82ycrc1CnE/bU1At+v0t0gw3SOq6OU11OiLDcEp3OXAWqqszllqkNy7O9n+BQ4mNet1o97ok91MiJSrSQChyL3BDNBxc6Z2XSDQcTMvBbwltnJYiX3EfyqBM+zIMw0jG2v33hwMOcIG/+/a1wN9+UC09Dp1Tev62SKouwo1YE9OH4Uax2eQ64EJUn/BGoTE+whk3+aa0t2RPtmjoPdOTzzAMoxISC/x97LHwxBNw5pn5lsjIHTtTdItOjBU450O+KVGRRsPBQ5KdG4ZhVGmOOgp23925DTzjDPN2VHVZAOyC24YTTzfgp3Qa8vWGBEKRmwOhSL0k6XUDocjN6XRoGIZRoYkF/v7mG4sKU7UZBdyLSGw9thUi5+DscdIy3fb7U+sWoEGS9HoUj9hiGIZRuTntNGjTxu0pzaqfdKPC4JwuPAe8jtsh8hbO+OhhVB9Ipym/ilRwEV8S6QiUKX6bYRhGhaVmTbj2WvjwQ3j33XxLY+QK1QG4nSL7AV2A7VC9CZFW6TRTqklaIBRZjlOgCswJhCLxyrQQqIPT4IZhGFWLPn1g8GC3VtqtW76lMXKF6irAOeUXaY7I7cB5uO2Wvkhl2305bjQ6BhiA83kYYx0QjYaDyayeDMMwKjf16jnHDAMGwBdfQIcOeRbIyAoijXH+CHrhnDuEgftwMa6vB2bhFKlvSlWk0XBwHEAgFPkZ+CAaDq5PW2jDMIzKyqWXuhFpOOy2wxhVgdtwlrnjgCNwQcMPw62THonq2+k26He38ddAw0AouQVbNBy0dVLDMKoejRvDJZfA8OEu8HfbtvmWyMicINAH1amIPAjMBn5CtV9ZG/RrbPQHzm1fSYdhGEbVpF8/Z3x05535lsTIDtvjpm/xotSsAUZn0qDfEWmiQ4aaOIvdS3B+EQ3DMKomLVrAuec6/7uDBuVbGiNzCnBrozE24gKGlxm/0V+SzRlPDYQic4ALsAgshmFUZa69FkaPhrvvNmf2lR8BHkNkrXddBxiNSFFlqnqM3wYz9cj8BW7R1jAMo+ryj3/AKafAQw9R0Lt3vqUxMmNcwvVjmTZYZkUaCEUa4IJ6/5KpEIZhGBWeUAieeIKG48fD3nvnWxqjrKj2yXaTvhRpnGOGGIJzD7gSsJ9nhmFUffbeG448kkaPPAJDhrh9poaB/xHp5QnXm3DWuh9Hw8Fl2RXJMAyjgnLDDRR26wZjxsDliV+LRnXFr7FR4pyyYRhG9ePgg1nTuTN17rwTLrrIbYsxqj2+10gDoUgd4AygvZc0C3g8Gg6uzoVghmEYFZG/Lr6YOhdc4DwdnXVWvsUxKgB+10j3ASbhnPh+7SWfBwwNhCLBaDg4I0fyGYZhFGPqrEW+yy5evJzZq/yVjy/bs32zpGVWH3oo7LGHcxvYu7cF/jZ8ezYaBbwH7BANB7tFw8FuQCvgHS/PMAyjeiDiLHhnzYKXX863NEYFwK8i3R0YFA0HV8YSvPMhXp5hGEb14dRTIRCwwN8G4F+RfofzT5hIC+CH7IljGIZRCahRw3k7+vhjeDvtYCFGFaPENdJAKLJ13OVA4N5AKDIE+MhL6+Klh3InnmEYRgUlPvB3jx75lsbII6UZG/1BcScME+PSxPv7IlCYfdEMwzAqMHXrwtVXww03wIwZsM8++ZbIyBOlKdLEiC9Gjom3REzH0tAP2W7PMAxcrNJhw+D22+HJJ/MtjZEnSlSkJUR8MQzDMGJstZWLBnPHHfDjj7DzzvmWyMgDpa2R7gN8EQ0HN3nnJWL7SA3DqLZcdZULr3bHHS7UmlHtKM1qdzqwbdz5p97fxOPTXApoGIZRoWneHM47D8aNg/nz8y1N/hC5FJGfEVmDyGeIHByXdxciSxH5BZHeCfWORuQ9RCSxycpCaYp0R5xj+tj5Tt7fxGOnXApoGIZR4bn2WtiwAdq1g4ICdjjoIJgwIfv9TJjg9q8WFLi/ueijLIicCtwD3AZ0BD4AXkGkNSJH49zL9gKuA/6LyLZevYbA3UBftPJuyC1tjXQuQCAUqQlcBjwQSzMMw6jqlOSGMJnhXrNJr7B7QQEFy5cDUGP+fDaefz4/zfiOJd0PK1J2ydKlLNx6HgAHtN0W30yaBLfcAmvWuOu5c6FvX3ee/2Dj1wBjUY3NbV+ByBHAJcAyYBqqbhZTZARuEPYHTvE+huqsPMicNVL62o2Gg+sDocilwIPlII9hGEalo+2IYRRs3FgkrXDtWna561a469bcdbxqFQwYkF9FKlIL6AQMT8h5DTgQpyz7ItIEN4NZF5iNSBfc7pBKv2/Ib/SXKcChwJgcypJTtt56a6ZNm1amuuvWrSMajWZVnmSsXrNh83mdDetZvSZ721Uyba8s9dOp46dsqjKZ5lcGyvse7D30V6bOb8nXRhX47IaBRdI2btxAYaH76q1b07/D+/ZDhpBsEVHnzePtMn63+aSGiEyPux6lqvE+1rfF+RJIfHCLgJ6oTkHkMZw9zWrgHGAFzk/7xUAfRPoBq4ArUP0gN7eRO/wq0jeA2wKhyF7AZ8DK+MxoOPhctgUrDW+EfC3OReE3QL9oOPhuaXWWLl1KjzJ6H4lGowQCgTLVTYf4qaTlixfTqGnTrLWdaXtlqZ9OHT9lU5XJNL8yUN73YO+hvzJrmrek7sJfi5Vd02IH/jyzaADwxYsX09Sr37mECDNJGTfOTecmIK1bl/m7zScbVLVzRi2oDgIGbb4WGYBbR/0L57O9A7An8BQiO6G6LqP+yhm/ivR+7++VSfKUcvRsFAhFYoval+Ii0lwKvBIIRdpHw8F5uejzg6g5MzAMo2Rm97uB9rf0p3DNlvDMG+vUZXa/G7LXydChbk101aotafXqufT88gewEUj8VdAM+K1YaZFdgPNxRknnAO+guhBYiEhtYFe2hOusFPhSpNFwsCIF3LsGGBsNBzcvagdCkdiidhbfWsMwDH8sOupEwK2V1vltPquaNmfONQM3p5dEOnFV6diTZrfcubmPNc1bMrvfDSzq2BNStFNSbNWsoLoOkc+Aw4Cn43IOA54tUtZtcRkF9Ef1L0QKgJpxeTWphC5n/Qb2Pht4MhoOrk1IrwWcFg0HH82FcEnkSLWoXQQR6Qt4Zm2oiKxOLOOTQtwvrvIk231m2l5Z6qdTx0/ZVGVS5dcANpSSXxko73fR3sOylFm0sJDrL9vI9Zcly83Oe7jwV7j+MkroI5vU9VHmLmA8Ip8A7+PWPrcHHk4odz6wDNXYcuB7wBBEDgL2BtYD32dF6nLE79TuI8CrwOKE9IZeXrkoUlItaifgLYhnHHhcREapat/UJbNHtvvMtL2y1E+njp+yqcr4yJ+e8VpPninvd9Hew+yXqQrvYTFUn0RkG1xEsBbATOD/UN2yqCvSzMvvGldvOiLDgOeB5cBZqJZ1wJM3/CpSoWgkmBitcYvFVZ2Xq0CfmbZXlvrp1PFTNlWZfHxO5U1536O9h7krU7VQfZDStkmqLgICSdKHAcNyJVZ5UKoiDYQiX+MUqAJvB0KR+OmIQqANMDl34hUjvUXtLKGq5f5Pke0+M22vLPXTqeOnbKoy+ficypvyvkd7D3NXxqg6pBqRPuP93QOI4Pb+xFgHRElcTM4h0XBwXSAU8beobRjFyXia3zCygL2HVQzx494wEIqcAzyRaGyUD7ztL+Nx215ii9rnA7ubC0PDMAyjvPG7RjoZaITnxD4QiuwJnAp8Ew0HH8+RbEmJhoNPBkKRYovapkQNwzCMfOB3RPoWMD4aDo4JhCLbAj8CC4AdgCHRcPA/uRXTMAzDMComfh0t7AV85J2fBMyOhoO7A2cDF+VCMMMwDMOoDPhVpHXZYmjUE3jJO58BtMq2UIZhGIZRWfCrSH8ETgiEIq1wwVlf89KbAX/mQC7DKBdEpJWITBORWSLylYicnG+ZjOqHiDQWkeki8oWIzBSRC/Mtk+Efv8ZGg4HHgf8Ab0TDwY+99MOBz3MhmGGUExuAfqr6hYg0Bz4TkcmqujJVRcPIIsuBbqq6SkTqAzNF5DlVXZJvwYzU+BqRemHSWgOdgSPisqbinMgbRqVEVReq6hfe+W84px9b51Uoo9qhqhtVNRbWpTbOm1yy8KNGBcSX1a5hVFREpBvQHxfMYHugj6qOTShTLH6tqhaLXysinYBxqrpHruU2qhbZeA9FpDHwNrAzcK2qPlAuwhsZU+LUbiAUuRe4IRoOrvTOSyQaDiaLU2oY5UED3F7iR0kSPEFEksavFZH2qjovrtzWXn1bmzLKQsbvoar+Cewtzrn7cyLyjDr/tEYFp7Q10j2JxYlz5yVhQ1ojb6jqZDx/zyIyNkmRa4Cxqro5fq2IFIlfKy6Y8AtAWFU/yLXMRtUjG+9hXFuLRORL4GC2uGk1KjAlKtJoOHhIsnPDqCyISMr4teKCCY8F3lTV8eUqoFEt8PkeNgNWqepyEdkK6AY8VK6CGmXG7/YXw6iMlBa/trl33hXn7vI4b+vBFyJS2gyMYaSLn/ewDfCuNxJ9F7hPVb8uPxGNTEi5/SUQitQFrgNOBHbCTeXOwUVf+U80HKx0QVgNI4aqvof9oDTyjKp+AnTItxxG2Sj1CyQQitQA3gRuBH4G7gMeAOYCNwNTvTKGURHJS/xaw0jA3sMqTiol2BdoC+wTDQe/ic8IhCJ7AG/hrBxtLt+ocKjqOhGx+LVGXrH3sOqTakrrJGBoohIFiIaDM4FhgLlUM/KGiDQQkQ4i0gH3Prf2rlt7Re4CzhWRC0SknYjcg9vn93CeRDaqIPYeVm9SKdLdcVO7JTEVsM3rRj7pjHNT+TkuuMJg73wIgKo+CfTDxa/9AjgI+D9Vtfi1Rjax97Aak2pqtwleMO8S+B1onDVpDCNNVHUaKVypqeqDwIPlIpBRLbH3sHqTakRaiHPqXRKbvDKGYRiGUS1JNSIV4LFAKLK2hPzaWZbHMAzDMCoVqRTpOB9tFPMraRiGYRjVBYv+YhiGYRgZYB5dDMMwDCMDTJEahmEYRgaYIjUMwzCMDDBFWsUQkbEiMilFmYCIqIh0TqPdQSIyM3MJqxYiMk1E7i+HflRETsp1PxUVEYl6z0BFpHnqGpvr9fDqbJtFWXx/5iJybpzcOX9PjPxgirQC40cpJuEq4My4NpL90/8CtMB5WMkKnnL+n4jMEZHV3t9hIlI3oYzGHStE5HsR+a+I7JUtWcqZE0gIzJwJpXzmLYCXs9VPJWUI7jksTqPOB16dJTmRKDVPev1/mKf+jXLAIrdUMVT1Lx9lNpL9qBO74ZxzXAL8CLQDRgHb4IIfxHME8CXOlVo74GLgMxE5S1Wf8NuhiNRS1XVZkL3Mbavq0lz0n6SfChklpKTnJCI1VXV9Gdorrd7ydJ+DJ1venp2qrgZWi0hO3lOjYmAj0kpEbLQiIleJyHwRWSYij4hIvcQysXOgO3BZ3CgwkDi1KyKF3mjyZ280+aOIXCcivt8PVX1VVc9V1SmqOkdVI8BQXBzbRJao6m+q+rOqTlbVY3BRMR4Wkcal3H/Um2IeIyJ/AhO89ANF5G0RWeU9l4dEpFFcvWki8rCI3OM9s2Uicmf8/ZXS9gki8rWIrBWRX0RkgIhIQtv3x13XEpHbReRXT55PReTwhPvYTUReEpG/vFH5hyKyp4gMAs4BgnGfVw+vTpGpXa/8VO/zWup97lvF5ad8V0p4xu1FJCIiy0VksYg8LnFTqXHtXi8ivwK/xr1Pp4vImyKyGrhIRApE5Cbvua31nuOxcW0lrVeafAmyxqZtjxIXkH2NiHwmIp2SlNnWu/6fiHwj3kyJ9+6/K3GzACJytNfOGu9/YqiI1CpFjhNE5Ku4z+JtEUkMmWZUYUyRVj4OxgUK6AmcChyPm85NxlW4KaVHcNNLLXDTuokUAPOBU3AjxAG4GLR9MpS1EbDMZ9nhwFa4+yqNa4DvcE7CbxSRPYHXgJeAvXFTrR2AMQn1euPu8wDcl3VfnBPx0truhFPwzwF7AiHcNO7lpcj3CO7Hyxm4z2kc8LKI7A0gItsD7wGKC6O1Dy7Gb6H3DJ7CBYOIfV4fJHYgIvWBKcAKYD/cO3BgkntO511BRFoA7wAzvXZ7Ag2AFxN+VHUH9sLNLPwzLn0Yzpdse+AFr69rgetxz+954DlxEVIopV66DPf66AzMASaV8oPhSqCmVwfcu74zcB6A96NnAnA/LmjHebgoWLcla8z7kfEE7nNuB3QDxpfhHozKjKraUUEPYCwwKeH6F6AwLm00MLWUOtOA+xPaDeC+yDuX0nc4od1BwMw0ZG+DC2h8jZ9+gTpe3nWltBkFXk5IexT4X0JaB6+tpnHP4Ac8ByRe2kDg1xRtTwDeTEgblFBv8/MF/oHzP906oc4LwIPe+VBgLlDLz2cel67ASd75hcBfQMO4/B5embZ+35UkfQwB3khIa+K1u19cu78DtZN8rv9KqDsfuDkhbRrwWGn1Svns+yekxe65d1xaA+BP4IKEMtvGldkXWOfd73rgyLi8d4CbEvo5DvejReLuIfaZ7+O13yaF/Jvr2FH1DhuRVj5mqVvjjLEAaJppoyJysYhMF5HfRWQFcDXQOlW9EtpqBrwKvA7c7bea9zeVq63pCdedgDO9KdIVnuzve3n/iCv3kXrfaB4fAi0lbgo4Sdvt4tqK8V6SejH28e5jVoI8wThZOgLvaWZru+2Ar1R1eVzaBzgl3j4uLd13pRPQLUH22AxG/LOcqarJ/G9vfn7e89me5M+vfUJa4nNPl82GPKq6Avg6SR/ElfkU94PmJmCUqr4Sl90JGJDwDCYC9YFk1sJf4mYQZorIsyJyiYhsl+H9GJUMMzaqfCQaYigZTtGLyKnACKA/7gv5b+Ay3FRgum01x8WwnQmclaC8SiP2xTcnRbmVCdcFwH9JrrDn++y7pLZLI9l9FXjp+1L8c1qdpixlJV6udN+VAiCCew8SWRR3XtJz8vv8Ep9dOs89Y7w17oOAjcA/RETi3tMCXCzRp5NULRZSUlU3ikgvoAvQCzgfGCYi3VX1y5zcgFHhMEVa9VlH6lB3BwEfq2q80cw/SimfFG+N7S3gG+B0VS0tBF8i/XHTlVPT7HYGsLuqzk5Rbv+EL8wuwAJV/buUOt8CXRPSDsJN7S5PUv5z3Ii0uaq+VUKbn+NG0CVZBfv5vL4FzhORhnFyHIhTAt+mqFsaM3Dr5HO1DBa38ajq3yKyAPf83ojLOgiYlUnbSeiC9wPMWz/eg9KDaVyDmz3oBkwGrgDu9fJmALv5eJ82471THwIfisgQ3Pt/Km60alQDbGq36hMF9vMsJLeV5Ja4PwD7iMiRIrKziNyEMyjxjWdE8zZuq0E/YFsRae4diYphGy99R6/Pl3AGHRerj+07Cdzu3d/DItJRRNp6VpwjE8ptD4wQkV3FWb9eS+pp5/8A3cVZ8+4iIr2BfwF3JCusqj/g1lXHishJIrKTiHQWkf4icoJX7EHcOt5TIrKvJ+/pcQY4UWAPT85tRaRmkq4mAKuAR8VZ73YDRgLPpaMAkvAAzuDrSRHZ35O/p4iMEpGGZWjvTqC/d3+7eErmYLYY+mSLgSJymIjsjjO4Woebji2GZ/Q1FLhQVT8ALgVu9+qCWzc9Q0SGiMge4iysTxKRpJ+5iHQRkYHeZ9kaOAZoRfZ/LBgVGFOkVZ/huC+WWbipqWTrniNx1qITgU9xRiD/SbOfXjjrx+7APGBh3NEqoeyrXvo3wD2eXJ01jT2kMVT1K9zIIoBT5F/irEAXJRSdgBvpfYwzuvkfKRSpqs4ATsZt4ZmJM8AK4yw6S6IPznL3DpwF8CRPvrlem/O961q40fvnuBFRbPQ+GjeqnI57LokjYlR1FXA4zir6E+BF3IjovNLuJxWqGhtBbsJ9Rt/glOta70iXe3HK9A7c8zseODEHU54h3Ps6A/cOHqWqxaaLRaQO7j2YqKrPAajqROAZYKKI1FbVKbg17UNwz/YTr/15JfT9F+6ZTcLtn/4P8G9VfSx7t2dUdCyMmlHlEZFpOAOZ0ratlLXtD4G3VTWU7baNLYhIFGf1OjwurQfux8h2qvpHfiTzRy7fQSP/2IjUMMqAiNQW59Bid9xoy8g9Qz1L2oyt1MsLEentWf4enG9ZjNxhxkaGUTaOxBm0vITzp2rklu44RwqQP7+5ZeEl3HICuP2tRhXEpnYNwzAMIwNsatcwDMMwMsAUqWEYhmFkgClSwzAMw8gAU6SGYRiGkQGmSA3DMAwjA/4fm+awAx2exVMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x216 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "mpl.rc('font', size=14)\n",
    "print(f'Recall of oracle @1m: {np.mean(err_init<1)*100:.3}%')\n",
    "\n",
    "fig, ax2 = plt.subplots(figsize=(6, 3))\n",
    "ax2.set_xlabel('Initial 2D reprojection error [pixels]')\n",
    "plt.grid(axis='x', which='both', alpha=0.5)\n",
    "\n",
    "color = 'tab:blue'\n",
    "ax2.hist(err2d_init, bins=bins, alpha=0.3, color=color, zorder=5);\n",
    "ax2.set_ylabel('Distribution of initial errors', color=color)\n",
    "ax2.tick_params(axis='y', labelcolor=color)\n",
    "ax2.yaxis.set_label_coords(-0.11,0.42)\n",
    "\n",
    "ax1 = ax2.twinx()\n",
    "color = 'red'\n",
    "ax1.plot((bins[1:]+bins[:-1])/2, prob, 'o-', c=color, zorder=20)\n",
    "ax1.set_ylabel('Recall ($\\Delta t < $ 1m)', color=color)\n",
    "ax1.tick_params(axis='y', labelcolor=color)\n",
    "ax1.set_xscale('log')\n",
    "ax1.set_ylim([-.03, 1.03])\n",
    "ax1.grid()\n",
    "ax1.yaxis.set_major_formatter(mpl.ticker.PercentFormatter(1.0))\n",
    "\n",
    "# plt.savefig('convergence_stats_Aachen.pdf', bbox_inches='tight', pad_inches=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Cumulative recall"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Prediction: [70.51142546 75.40805223 81.93688792]\n",
      "Initial:    [ 0.10881393  1.19695321 47.8781284 ]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ0AAAEnCAYAAAByjp6xAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA9nklEQVR4nO3dd3wUdf7H8dcnCUkIofdeBQFBQESagr17inr2A69g/+nZu9x51vPseoL9bIeKvaBwGhELvfcWQk0CgZBC6n5+f8ygSwzJJuzu7G4+z8djHrs7Mzv7TobNh5n5zvcrqooxxhgTDnFeBzDGGFN3WNExxhgTNlZ0jDHGhI0VHWOMMWFjRccYY0zYWNExxhgTNlZ0jIkAIvKaiBR5ncOYULOiY8wBiMgfRURFZJXXWYyJFVZ0jDmwS4F0oKeIHOlxFmNighUdYyohIh2AUcDtwGacAmSMOUhWdIyp3MVAIfAJMBm4QETi/VcQkXEiMl1EtotIsYisEZE7ROQ33ysRGSwin4pIjogUishSEbmjkvXai8hHIpIvItki8lglnysicp2ILBGRIhHJEpGXRaRFhfXSRWSqiIwUkdnuuutF5A/B+AUZUxtifa8Z81sishhYrKqXisgRwFzgVFWd6rfOHGAFsBAoAo4HxgCPqOrtfusdD3wOZAGvAluBXsAwVR3mrvMaTqFbBcwG5gAnAOcCV6vqv/229wLwJ+B1d72OwHVABnCkqha566UDxUBj4GX3c/8IDAT6qeqyoPyyjKkJVbXJJpv8JqA/oMAZfvNWA29WWC+lkvdOAvKBJPd1HLAO2AQ0q7Cu+D1/zf3MeyusMx+Y6/d6uLveHyqsN9KdP95vXro77xi/eS1xCuRjXv+ebaqbk51eM+a3LgVygK/85r0DnC0iDfbNUNVCABGJF5Gm7umt74AGwKHuaoOAbsBTqprj/yGqWtlphhcrvP7eff8+v8cpalNFpMW+CVgJZALHVnj/alWd4feZ2ThHU90wxgMJXgcwJpK412MuwikenUVk36LZOMXkbOAtd92RwIPAUUBihU01dh+7u49LA/j4UlXdVmHeLqCp3+ueQCpOgalMqwqvMypZp+I2jQkbKzrG7G800MGdzqlk+aXAWyLSDZiOc9rtrzh/3ItwjmweoXaNdHwBrBMH7AQuPMDyXRVelx9gPTnAfGNCyoqOMfu7FNgBXFXJspOBcSLSCjgLSALOVNWN+1YQka4V3rPOfTwMmMrBWwecCPysqvlB2J4xYWXXdIxxiUgyTmuxL1T1/YoT8C+c/6hdyK9HEOL3/iTg2gqbnQ+sB64XkWYVPq82RxuTcb6391aSP15E7LSZiWh2pGPMr84CGuHcm/MbqrpSRNbgHA1dCpQAn4nIRJyjnsuocIpMVX0iciVOk+mFIvIKTtPlHsAIdwqYqs4QkeeAW0SkP05jh2J3e+fhFKPXarJNY8LJio4xv9pXSL6uYp2PgZvd52fjNCR4FOeU3H+AtIrvV9VpIjIKuA+4EYjHOfp5szYhVfVaEZkPXAk8AJThXFN6F/imNts0Jlzs5lBjjDFhY9d0jDHGhE1Yi46IHCMin4jIFrfL+HEVlouITBCRrSKyV0TSRKRvhXWaisgbIpLrTm+ISJNw/hzGGGNqJ9xHOqk4N8ldD+ytZPmtwE04/UgdidNX1TQRaei3zts490Kc4k6DgDdCmNkYY0yQeHZNR0TygWtV9TX3teC06nlWVR9w59XHKTw3q+pEEekNLAdGquoP7jojcboKOVRVbbAtY4yJYJF0Tacr0Aa/lj+quheYgdPJIcAwnH6nfvR73w9Agd86xhhjIlQkNZlu4z5W7FMqE2jvt062f0eJqqoikuX3/v2IyHhgPEBycvIRnTp1CmpoEx4+n4+4uEj6P5KpCdt/0Wv16tU7VLVlsLYXSUUnJFR1Ek538/Tq1UtXrbIzcNEoLS2N0aNHex3D1JLtv+glIhurXytwkfRfj+3uY+sK81v7LdsOtPTvPsR93spvHWOMMREqkorOBpzCceK+GW5fWEfz6zWcn3BawA3ze98wnC7n/a/zGGOMiUBhPb0mIqk4fUSBU/A6icgAIEdVM0TkSeBOEVmJ02X83TgNB94GUNUVIjIVmOheqwGYCHxmLdeMMSbyhftIZzCwwJ3qA39zn//dXf4o8ATwHM6Y9G2Bk1Q1z28bFwOLcDo6/Mp9flk4whtjjDk4YT3SUdU0qhg8ym2VNsGdDrTOLpyOGY0xxkSZSLqmY4wxJsZZ0THGGBM2VnSMMcaEjRUdY4wxYWNFxxhjTNhY0THGGBM2VnSMMcaEjRUdY4wxYWNFxxhjTNhY0THGGBM2VnSMMcaEjRUdY4wxYWNFxxhjTNhY0THGGBM2VnSMMcaEjRUdY4wxYWNFxxhjTNhY0THGGBM2VnSMMcaEjRUdY4wxYWNFxxhjTNhY0THGGFOp5Vv3BH2bCUHfojHGmKiUtaeIz5dsY/HmXOak57B5196gf4YVHWOMqcNKy318szKLJ6atZuX2PADaNk6mV5uGjBnUgZseCe7nWdExxpg6avOuQs77909s31NE+yb1ueu03ow8pAW92zb6ZZ2bgvyZVnSMMaaOUVWy84r5v3cWkFNYwsNj+nFqv7Y0rl8v5J9tRccYY2LM+ux8svOK2VtaTlGpj+KycvaWlDN34y7WZeezansehSXlADx4Tj8uHNIpbNms6BhjTAzw+ZS3Zm3krVkZv1ybOZBxw7vQpXkKg7s047D2jcOU0HHAoiMiNW0rp8Dhqpp+UImMMcYEJCuviHnpu/h5/U6mLttO5p5i+rZrxF9P6MngLk1JrhdPcr04kuvFU79ePMn14mmUnEBCvHd3y1R1pJMK3ADkBrAdAZ7H7vsxxpiQytxTxD+/WsXsDTlk5BQCkJQQx6ieLTm9f1tO79fW06JSnepOr/1XVbMC2ZCIPBOEPMYYY/zsKijhs8VbmbFmB2uz8tmwowCAE3q34rKhnTmiS1MOa9eYxITILTT+Dlh0VLVGP4GqNjz4OMYYU3cVlZYzY3U267ILWJedz0/rdrJlt3ODZufmKfRt14gzD2/H4M5NOaZnS4/T1k5ENSQQkXhgAnAp0BbYBrwFTFDVMncdAe4DxgNNgVnANaq6zIvMxhhTU+k7Cli2dQ95RaXkFZWRkVPInPSc/RoAtG6URL/2jRk3vAt92zViWPfmOH/+oluNio6IpAL3AMcB8cBMnIKQE6Q8twHXAGOBJUB/4HWgGLjfXedWnPuVxgGrgHuBaSLSS1WrbrJhjDEeUVWmr8jisa9WsSpz/z9VIjCwY5NfCszJh7WhUXLo75nxQk2PdF7AaaV2H5AIXAW8AZwepDzDgU9V9VP3dbqIfAIcBb8c5dwAPKyqU9x5Y4Es4GJgYpByGGPMQdlbUs6qzDyWb93D8m25/LhuJ+uzC0iIE64e3Z0jOjeld9tGNExOoEFiAnFx0X8UE4gqi46I/FlVX/KbNQw4RFV97vLlwJwg5pkJXC0ih6rqShHpg3NU9ZC7vCvQBvh63xtUda+IzMApWFZ0jDFhlbu3lOVb97CrsIRdhSUs3ZLL3HTnJkyfOus0TE6gT9tGnDuoAxcc2ZEWqUnehvZQdUc6p4jIpcBfVHUN8DPwsoi8C9QDrgR+CGKeR4CGwHIRKXfzPaCqz7vL27iPmRXelwm0r2yDIjIe5/oPLVu2JC0tLYhxTbjk5+fbvotisbb/dhX5mLqhlOU5Prbl+yjTX5clxkGf5vGc2a0enRrF0alhHC3qCyLFwGaWzt3sWe5IUGXRUdXzRORs4EsReQW4DrgD+AfOPTkzcS78B8sFwB9wTpUtAwYAT4nIBlV9uTYbVNVJwCSAXr166ejRo4OT1IRVWloatu+iVzTvv+KycrbnFpFT4BzJzFqfw8szN1DmU47p2ZJT2zRkePfmtGqYTNMG9WjWIJGkhHivY0esaq/pqOpHIvIN8DDwLXCFqt4Sojz/BB5T1f+6r5eISGecQvcysN2d3xrI8Htfa79lxhhTI6pK7t5SMnIK2bizkOXb9pC+o4BtuUWszvy1n7J9erZO5ZFz+zOwU1OPEofJz/8O+iYDakigqntwrrUMByaJSBpwh6oWBDlPClBeYV45v/Z0sAGnuJyIey1JRJKBo4FQFUJjTJQo9yml5T7KfUpZuVLq81FWrpT5P/qU0jJlW+5eflq/k/kbd7F+RwF5RWW/bCchTujcPIW2jetz7qAO9O/QmOapiTRNSaRjs5S6cU2mKBem3hH0zVbXkKAT8BjQG1gM3AwcgXPksVBE/qqqnwUxz6fA7SKyAef02kDgRuA/AKqqIvIkcKeIrARWA3cD+cDbQcxhjIkQ+cVlbNu9l48WbmH+xt0Ulpazt6SMwpJyikqd3pNLyp1iolr99vyJwFFdm3H2gPZ0bp5Cp2YpdGqeQpfmDUiuV8dPkc19BaexcnBVd6TzH5wji1uAk4GJqnoW8He3McFEERmrqucHKc91OPfjPA+0wrk59EXg737rPArUB57j15tDT7J7dIyJbvnFZWzKKWRTTiEZOYVs3rWXmWudrl/26d+hMU1SEmnbKJmUxHjqJzodWSYmxJEQJyTEx5EQL87zuH3P/ebFu+vFCa0bJdO5eQpNUhI9/Kkj1O4MmPkkdBoOTA3qpqsrOoNxeo5eJyJf4ZzeAkBVVwKj3NZhQeEWjhvc6UDrKE7jhQnB+lxjTO0Ul5WTtaeYUvdIw/8UVlm5UubOX5hVRsmy7c58n5JfVMamXYW/FJlNu/aSU1Cy37ZTkxLo2CyFk/q05vT+beneMjXs3fDXScV58OppULQbTriPcBedeThHNa8DJ+D0ErAft3WYMaYOKPcpuwtL+GTRVj5fvI3FW3IpKfMF9ub58/Z7mRAntG9an07NUji5XWM6NUuhY7P6dGzqnOZqklIvJrp9iTrfPgS5m+Dcl6HT0KBvvrqi8wfgX8ATwELgiqAnMMaEVVm5j50FJWTtKWZPUSl7S8r3u06y71rJvuf75q/KzGNTTuEvNzx2a9GAPwztTM/WDZ3TW36nteLjhXp+p7UWL1zAkCMH/3K6KyUxntaNkomvI3fhR4010+Dn56DP76DfeSH5iOru09kIhOaTjTEhk7GzkHkZOWzdXcTO/BLyi52OJXcXljIvY1e1RycikFIvnvqJCdRPjCOlXgKHtErlrMPb0bxBIi0bJnN871YBX2zPT4+3U2ORrrQIvrgZkhrDWaEbqSaiepk2xtSMqlJYUs6uwhLmpu9ixfY9zFyzg2Vbfx34t0FiPA2T65GanEDD5ATGDGzPYe0b06phEo3r1yPFLSz1ExPcQhNPUkKcndqqS5Z9CF/eBvmZcN6rkBy6/yBUNVx1IdBZVbMD2ZCIZAFDbLhqY4Irt7CUr5ZvZ3tuEVl5RSzctJusPcUUFJdRWFq+XzPhhDhhUKemXDmqOyf2aUWfto2pn1jHm/6aqqnC9L9BQhJcMgUOOSGkH1fVkU4ycIaIBDJcNUADbLhqY4Jie24RczfmsDYrn9d/TGdXYSkATVPq0bN1Q47v3ZiUxAQaJMbTICmBBkkJDOjY5JfrK8YELH0m7NoApzwc8oID1Z9eq1V/Z8aYmikr97Fyex6rM/P4fs0OPlywBXCurXRqlsJDY/px3KGtraCY4PvpWUhuAodfFJaPC9pw1caYmikuK2du+i4mz9lE2qos9rjdsCQmxPH7wR0YM6gDAzo2sTvjTehsmgOrp8LRN0P9JmH5SGtIYIwHisvK+f3En1m0aTcNkxM4pW8bRh7Sgr7u/Sp2RGPCYtHbUC8FRt4Qto+0omNMGGXnFXPjuwv5Ye0OfAqn9WvDQ+f0p3FKbA5NbCJYUS7Mf8O5JyepYdg+1oqOMWGgqrz0/QYen7aaknIfY4d3YXj3Fhx/aKs6M0yxiTBpj4CvFI78U1g/1oqOMSG2fOseJny6jNkbcjju0FbcdXpvurdM9TqWqcvys2Heq3DISdB5eFg/2oqOMSGgqvy0fiev/pDO9BWZJCXEccvJvbh6dHe76dJ4b+bjUFYMJz8Y9o+2omNMEJX7lNd+TOfdOZtYlZlH05R6XD26O2OHd6FVw2Sv4xnjSJ8JXY+BFoeE/aOr6pEgjwBH8FHVRkFLZEyU8fmco5q56bv4z0/p7CwooXvLBjx6bn/OGtDOmjybyLJ5HmStgKFXefLxVR3pXBu2FMZEsbs+Wso7szMAaNc4mT+N7Mrdp/e202gm8qjCJ9c6zaQH/9GTCFXdHPp6OIMYE40mz8ngndkZnNC7NU9dOIAGSXbG2kSwzXMhazmc8SQ06+pJBPuGGFNDuYWlrMnKY8r8LbwzO4OBnZrwxAWHW8ExkW/pFEhIhsPO9SyCXdMxJgCqykcLtzAnfRcfzt/C3tJyAI4/tBX/PP9wGibbzZ0mCmyY4YwGmuzdn2y7pmNMNT5bvJVnv1nLyu15NEiM59hDW3L+4I50aFKf7i1T7eZOEx02/ghZy6DffZ7GsGs6xhyAz6c8MX01z3yzlg5N6/O3s/ryh2GdrYGAiU6LJzuPgy/3NIadhDamgpyCEtZl5/PJwq288fNGTu/flkfO7U+qXbMx0WrbYpj/H+hxItRv6mmUgL5FIpII3AVcBHQC9juBrap2I4KJCW/8vJH7P11OSbkPgN8P7sDDY/rbKTQT3f73N0hpDmMmeZ0k4COd+4ELgIeAJ4BbgC7AhcA9IUlmTBjtyC/miWmreWtWBt1aNOCeM/vQo2UqHZuleB3NmINTUuj0QHDE5ZDSzOs0ARed3wNXqupUEXkM+FhV14nICuBEYGLIEhoTQmsy87hh8kKWbd0DwLjhXbjr9N7Ui7fxbEwM8Pngx2egrAh6HO91GiDwotMaWO4+zweauM+nAo8EOZMxYXP3R0tZuT2Pu07rzaDOTTmis7fnu40Jqh+ehLQHodNw6Has12mAwItOBtDOfVwLnAzMA4YBe0MTzZjQyd1byms/pDNrQw73nNGHP4305u5sY0Jm+xL439+hVV+4/AuIkFaXgRadD4HjgZ+Bp4B3ROQvQHvgnyHKZkxIrMvO58JJP5OdV8zw7s255KhOXkcyJvjmvwESB5d9GDEFBwIsOqp6h9/z90VkEzACWK2qn4UqnDHBtCmnkOfT1jJl3hYaJifw1p+PYmi35sRbyzQTa2Y8BrMnOk2kG7b2Os1+anXjgarOAmYFOYsxIbNkcy7jXp1NTmEJFw/pxJWjulvLNBObVGHuq9Cyd0Q0ka4o0Pt0HgA2qeoLFeZfCbRXVWs2bSJSUWk5b83K4P7PltMoOYFPrx3JYe0bex3LmNDJWgF7NsOZT0dEE+mKAm0XehmwoJL584A/BC+OMcGjqvzlP3O5/7PldGxWn/euHG4Fx8S2olz45DqQeOh5stdpKhXo6bVWQHYl83fiNKc2JqLsKijhlvcX8/2aHZzStw2PX3A4KYnWjY2JcXNehi1z4dR/QsM2XqepVKBHOhnA0ZXMPwbYHLw4ICJtReR1EckWkSIRWS4io/yWi4hMEJGtIrJXRNJEpG8wM5jotj47n0temsW3q7K4enR3nr14oBUcE/uyVsCsidCsGxw13us0BxToN3Ei8ITbB9s37rzjcbrFCdrNoSLSBPgBmAmcjnN01Q3I8lvtVuAmYBywCrgXmCYivVQ1L1hZTHRam5XHhZN+prCknJfHDmZ0r1ZeRzImPGa/CEW74ZL3vE5SpUCbTP9LRFoATwOJ7uwS4ClVfTSIeW4Ftqmq/3WiDfueiNOn/A3Aw6o6xZ03FqcoXYx1x1NnzVq/k+fS1vHTuh2kJiXw0tjBDO/ewutYxoSHrxw2fAftBkLb/l6nqVLA5xxU9Q4R+QfQx521QlXzg5znbGCqiEwGjgW2Ai8Bz6mqAl2BNsDXfrn2isgMYDhWdOqcsnIfj361ipe+X09KYgKXDu3MVaO706phstfRjAmfBW/CzrUw4gavk1Srpie66+NcB1qoqsUhyNMNuBqnJ+uHgQHAM+6yZ3EKDkBmhfdl4vSO8BsiMh4YD9CyZUvS0tKCGtiER35+/m/23bZ8H88sLGJrvjKiXQIX906kQb1sls/L/qWjQBMZKtt/JjjEV8qwn+6hNKUDc3LbQ4T/ngO9T6ch8ApwLqDAIcB6EXkB2K6qE4KUJw6Y69cDwgIROQS4Bqfo1JiqTgImAfTq1UtHjx4djJwmzNLS0vDfdzPX7OD+N+aSlJDA0xcdxlmHt/MunKlWxf1ngmjuq1CaS+IFrzC6x3Fep6lWoK3XHsHp8HMQ+3fw+RlwThDzbIPf/Cd1Bc7AcQDb3ceKzbRb+y0zMS4rr4hLX55FalICX1x/tBUcU3eVFMC3D0LHodA9MoYuqE6gRecs4AZVXYhzpLPPCpxTYsHyA9CrwryewEb3+Qac4nLivoUikozTnPvHIOYwEcjnU+am5zDq0TQAJpzVl7aN63sbyhgvrfkaCrJg2DUR1alnVQK9ptMU50bQihoC5cGLwxPAjyJyFzAZGAj8H3AngKqqiDwJ3CkiK4HVwN04Y/y8HcQcJsJM3VDKNd98RUGJ88/tqQsHcFq/th6nMsZjC96E1NbQPfJPq+0TaNGZg3O086T7et/RzhUE8QhDVeeIyNnAgzjDYGe4j8/7rfYoToOG53CK4SzgJLtHJ3Z9tWw7/11VwqBOTThnYHsGdmpq3dkYszsD1v4PRt0KSalepwlYoEXnTuAr987/BOBG9/kQnF4JgkZVPwc+r2K5AhPcycS4hZt2c93bC2iTIrwzfihJCfFeRzImMix4y3kceKm3OWoooGs6qvojziihicA6nN4ItgLDVHV+6OKZuu7FGetpkBTPrUOSreAYs095GSx8C7ofC02iaxDCao90RKQe8CZwp6qODX0kY5xGA0/9bw1Tl23ntH5taZac63UkYyLH6qmQuwlOecjrJDVW7ZGOqpYCJ7F/qzVjQuqxr1fx1P/WcErfNvztLOvP1Zj9ZPwE8UnQ81Svk9RYoE2mPwDGhDKIMfts3b2Xl2ZuYMzA9jx78UCaNUis/k3G1CWZS6FZV4iPvt7TA02cAdwtIkcDc4EC/4Wq+niwg5m666npawC46eReSJTce2BM2Kz9H6xPg1G3e52kVgItOuOAXUB/d/KngBUdExTfrspi8txNjBnUnvZN7MZPY34j/XuIS4AR13udpFYCHdqga6iDGOPzKfd/upzEhDj+NNL+yRnzG/lZzg2hHY+CxBSv09RKoNd0jAm5135MZ/2OAq4//hD6trObP435jWn3QlEuHHeP10lqLfquQpmY4/MpD09dyaQZ6xnRozlXj+7udSRjIo8qrPsWep8JnYd5nabWrOgYTxUUl/HXyQv5enkmx/ZqyeO/H2CNB4ypzJppkL8depxY/boRzIqO8Yyqct8ny/h6eSbnDurAo+f1Jz7OCo4xlZo+AZIbQ9+zvU5yUOyajvHMze8t5v15m7l8RBf+9fvDreAYcyBL3oesZXDMLVAvult1Blx0RKS1iNwsIv8WkRbuvBEiYs2MTI39uHYHU+Zv5uKjOnHnab29jmNM5CraA5/fCM0PgSHjvU5z0AIqOiJyBLAKuAT4E9DIXXQi8EBooplYtS13L1e8MY/G9etx7xl9qBdvB9zGHNCKT5wWa2c/DwlJXqc5aIF+2x8DnlLVgUCx3/yvgBFBT2ViVlm5jyvemEdxuY9Hzu1Hcj3rOdqYKi2eDE27QocjvU4SFIEWnSOA1yuZvw1oHbw4JtbdNmUJizfncs/pvTnlMBv505gq/fgMbJgBh18YNcNRVyfQorMXZ5TOig4FsoIXx8SyTxZtZcr8zZzWrw2XDevidRxjItuW+fD13dD7LBj+f16nCZpAi87HwH0isu+EoopIF+ARYEoogpnYoap8uWQbt09ZTI9WqTw0pmL3fcaY31gzzXk87Z9R2+VNZQItOjcDzYBsIAWYCawFdgN3hySZiRkvfr+eq96aT4em9XnmooE0rl/P60jGRL7076FNP2jYxuskQRVoh597gJEichwwCKdYzVfV6aEMZ6JfWbmPZ75Zy4COTZhy1XC7F8eYQOzZ6gzUNvQqr5MEXUBFR0QGqOpCVf0G+CbEmUyMKCot59KXZpFXVMYVx3SzgmNMoJa8B74yGDTW6yRBF+jptfkislREbhORDiFNZGLGl0u3MXfjLi4a0pFTDoutUwTGhFT2KkhtAy0O8TpJ0AXa99qh/Hpj6AMiMhN4A3hfVXNDFc5Ep3kbc7j1/cWsyy6gVcMk7j2jr3XiaUygfOWwdjp0HOJ1kpAI6EhHVVer6n2q2hPnZtDFOD0RbBOR90IZ0EQPVWXSjHVc9OIsdhaU8NcTevLhNSOon2g3gBoTsI0/Qn4m9B3jdZKQqHEv06o6C5glIm8BLwCx+ZsxNbI9t4gbJi/g5/U5dGvZgP/+ZSitGiV7HcuY6LPgTUhMhZ4ne50kJGrU6ZWIdBWRu0VkBU6z6RzgzyFJZqKGqv5ScM4Z2J6p1x9jBceY2ijMgWUfwOEXQWIDr9OERKCt167BuaZzFLAUeAV4W1W3hDCbiRJvz87g5/U53H16b/58dDev4xgTvVZPhfIS6PM7r5OETKCn124D3gGuUNUlIcxjosjGnQXcMHkhCzJ2c1j7Rowd3sXrSMZEL1WY+SQ07QJdRnqdJmQCLTqdVVVDmsREld2FJdz47iJWbNvDvWf04bzBHWyIAmMOxrxXYccqOOOJmOncszIHLDoiMghYqKo+YGBVTV5VdX4IspkIpapc/tocFmTs5t4z+vDHkTaOnzEHJXMZfHErdBrmXM+JYVUd6cwF2uD0Ij0XUKCyyqOAtYmtIwqKy7j8Vafg3H16bys4xgTDR1c7DQd+91zUD0ddnaqKTlecDj73PTeGez5ayuz0HK47rgeXj7B/FsYctEWTYdtCOOURaN7d6zQhd8CT8Kq60e86jgIZ7rz9JndZSIjIHSKiIvKs3zwRkQkislVE9opImoj0DVUG86sNOwr4eNFWzjy8HTed1Mv6UjMmGBa+6TQeGHy510nCItArvxuAlhVnikhzd1nQichQYDxO7wf+bgVuAq4DjsQ5/TdNRBqGIof51T0fLSVO4JaTenkdxZjYsGcrZMyCrsdAQlL168eAQIuOUPkRTSpQFLw47oeJNAbeAv4I7PKbL8ANwMOqOkVVlwJjgYbAxcHOYX61cWcBM9fu4IpjutOpeewMKGWMp+a8DOXFMPhPXicJmyqbTIvI0+5TBR4SkUK/xfHAEGBhCHJNwulM9FsRuc9vflecxg1f75uhqntFZAYwHJgYgiwGeOabtQCc2s96izYmKMrLYOVnzkBt7QZ4nSZsqrtPp5/7KEBvoMRvWQkwH3gsmIFE5C9AD+DSShbv+4uXWWF+JtD+ANsbj3OajpYtW5KWlhacoHVIYany8YJChraNJ3v1AtJWhz9Dfn6+7bsoZvvvtzps+pge2StZ1udmsuvQ76bKoqOqxwKIyKvA9e4IoiEjIr2AB4GRqloajG2q6iScIyd69eqlo0ePDsZm65Tn09ZS6lvFHWOGcnjHJp5kSEtLw/Zd9LL9V0FZCTx9NaS2pu/5d8f0zaAVBTq0weWhLjiuYUALYJmIlIlIGTAKuNp9vtNdr3WF97UGtochX52zI7+YJ6et4aiuzTwrOMbEnKVTYM8WOP7eOlVwoAZDG4jIscBFQCcg0X+Zqh4XpDwf4dyI6u9VYA3OEdBqnOJyIjDHzZUMHA3cEqQMxlVcVs6t7y+mpNzHfWdaq3RjgiJvO0y9DVoeCofXvfZPgfYyPQ5n7JwPgdHAx0BPnAv7bwYrjKruBnZX+OwCIMdtqYaIPAncKSIrcYrQ3UA+8HawchjYlFPI2Fdnsz67gDtPO5Q+7Rp5HcmY2PDRVVCUC5d9BHF1r7/CQI90bgauVdWXRCQPuENV17s3beaHLl6lHgXqA88BTYFZwEmqmhfmHDErK6+ICyf9TF5RKZMuO4KT+lqLNWOCIns1rPsGjrkF2g/yOo0nAi063YDp7vNinPtzAJ4F0oDbgxvrV6o6usJrBSa4kwmBR6euYsvuvbx35TCO7NLM6zjGxI6NPziPA+reabV9Aj2224lzAybAFuAw93lznKMOEyNmrtnB+/M2c9GQTlZwjAkmVWco6kbtoWnd7bcw0COd74GTgCXAu8DTInIicDwwLUTZTJipKu/N20ScwD1n9PY6jjGxQxXeGwdb5sJZz9a5Fmv+Ai061wL7Br1/CCgDRuAUoH+EIJfxwP2freDjhVsZN7wLKYkBN2w0xlRn5zpY/hEMuBQGXOJ1Gk8F9JdFVXP8nvuAR0KWyHhiR34xb8/eyOn92nLPGX28jmNM7Cgvgyl/AgSOuqJOtljzV9XIoQGf0PcvSib6lPuUP742h6JSH1eM6mZDFhgTLD4fvHuZM17OWc9C2/5eJ/JcVUc6O6h+rJx9vU/byKFR7Otl21m8OZdxw7vQv0MTr+MYEzu+vBVWfQFH3wQDK+tOsu6pqugcG7YUxjNFpeXc/N4ierdtxM0n2zg5xgTNlnkw50U47Dw4tm71r1aVAxYdVf0unEGMNx6ftpqCknJuPqknqUnWeMCYoCjYCZ9eD3EJcMpDdf46jr9Au8Gp8vqOXdOJPqXlPm6fsoQp8zczZmB7ju3VyutIxsSOmY/D9iXOdZxU+275C/S/ttVd37FrOlFm8pxNTJm/mXHDu3Dnab2Js8YDxgTPlnnQYQgMuszrJBEn0KJT8fpOPWAgcBVOh5smSizdksuT09cwfUUmfds14u7Te5MQb4f+xgRNaRFsngODxnqdJCIFep9OZdd3povIeuDPWA/PUWFXQQm/n/gThSXlnH9EB+4+o48VHGOCbfNs8JVB16O9ThKRDvbK8ULgmCDkMCGWV1TKJS/NoqTMx4Qz+zBuRN3t+8mYkJrxGKS0gENO9jpJRKp10RGRVOAGYFPQ0piQUFWue2cBqzLzeGnsYGs0YEyo/PQ8bPgOjr4ZElO8ThORAm29lsf+DQkESAEKgLrdkVAUWJ2ZT9qqbG49pZcVHGNCZeUX8NWd0Gk4jL7D6zQRqyYdfvrzAdnALFXdFdxIJpiKSsu5bcpikuvFcdbh7byOY0xsKt3r9K+WkAznvADxds/bgQTakOD1UAcxwefzKTe9u4hFm3fz70uOoENTO9w3JiQ2zYbSQjhnEjTt7HWaiFajcuzeJNqKCoO/qeryYIYywfHWrI18vmQbt596KKccZkNOGxMSe7bBtHshqRH0PMnrNBEv0Gs6A4FXgX77ZuFc47EOPyPUF0u2cc/HyzisfSPGH93N6zjGxKbCHPjP72DHajhnItRv6nWiiBfokc4rOMNUXw9kUn3v08ZDe4pKeeDzFfRolcr7Vw633gaMCQVV+PT/YMcquOi/0OtUrxNFhUCLziHA+aq6NpRhTHC8/P0Gtuzey1t/PorkenYQakxIrPwcVnzq9DxgBSdggd6OPhPoHcogJjg27yrkmW/WkJqUwJCuAY/DZ4ypCVWYPRFS28Bpj3mdJqoEeqTzJ+AlEekGLAVK/Req6oxgBzM1V+5T/vHZCnwKb//lKOpZFzfGhMba/8GGGc44OQmJXqeJKjU5vTYQqKxfB2tIEAF8PuWW9xYxddl27jj1UBsB1JhQKc6Dr+6A1NYw7Gqv00SdQIvORGA68BDWkCAi/bhuJx8s2MJ1x/XgilHdvY5jTOya/x+ntdqFb0NiA6/TRJ1Ai04H4DRVXRfKMKZ2VJV35mSQmBDHNcf28DqOMbGrtAgWvAUN20Kv07xOE5UCPek/DTgilEFM7b36QzqfL97Gcb1aWWs1Y0Jp5hOQtQxG3ABityLURqBHOlOBf4lIf2AJv21I8EGwg5nA/G9FJn//bDnHHdqKZy8e6HUcY2JX1gr44UnoPBKGXul1mqgVaNF53n28s5Jl1pDAQ09/s5YOTevz3MWDbEA2Y0LF54N3LoTEVBgz0es0US3QDj/tr1kEUVXmZ+zm/XmbWLRpNxPO7EP9RKv7xoTMplmwKx3O/jc07uB1mqhm/W9HmZIyH1e8MZdvV2WTlBDHaf3acOGQTl7HMia2Lf/YGbag91leJ4l6gXb4eWNVy1X18eDEMdV5dOpKvnUHZPvDsC6kJtn/G4wJqW2LYP7r0HkEJKV6nSbqBfoX67oKr+sBbYG9QBZgRScMfly7g5dmbmDssM5cPdqaRhsTFl/e5gxbcPbz1a9rqhXQtRpV7Vph6gC0A2YANwUrjIjcISJzRGSPiGSLyKcicliFdUREJojIVhHZKyJpItI3WBki2bPfrqV5g0RuP9W6wTMm5Hw++PQGyPgJhvwFGtqYVMFQ6wYCqpoJ3AU8Grw4jMZpKTccOA4oA6a7g8ftcytOobsOOBLnSGuaiDQMYo6I887sDH5ct5Orj+1hjQaMCTWfD76+C+a9CgMvc+7LMUFxsBcE4oDWwQgCoKr79e0mIpcBucAI4FMREeAG4GFVneKuMxan8FyM011PzNm8q5C/f7qcod2acelQazRgTMh9fqNTcNoNhLOesRtBgyjQhgRjKs7CuaZzDfB9sEP5aYhT2Ha5r7sCbYCv962gqntFZAbO0VFMFp3nvl1HcVk5D43pT1KCHeUYE1K7Njr9q/W/EH73nBWcIAv0SOf9Cq8VyAa+IYjXdCrxFLAQ+Ml9ve+kamaF9TKB9pVtQETGA+MBWrZsSVpaWtBDhtLyneW8O7eIoW0T2Lh0Dhu9DuSR/Pz8qNt35lfRsv+6rXudDps/QSWeucmj2Pv9TK8jxZyIvTlURB4HRgIjVbW8tttR1UnAJIBevXrp6NGjgxMwDN6ZncGjU5fQvEEiD18ygk7NU7yO5Jm0tDSiad+Z/UX8/svPhvcvh03fQ8ej4JwXOKpZN69TxaSIvMlDRJ4ALgSOVdX1fou2u4+tgQy/+a39lsWE3L2lPPj5CoZ3b85zFw+iaQMbKMqYkPn6Lkj/Hk74Gwy92gZmC6Eqj2BE5FQRSReRRpUsa+wuOzGYgUTkKeAi4DhVXVlh8Qac4nKi3/rJwNHAj8HM4aWC4jJufX8RecVlXH/8IVZwjAmlJe/D4snQ73wYeYMVnBCr7kjnWuCfqrqn4gJVzRWRR3Bak00LRhgReQ64DDgb2CUi+67h5KtqvqqqiDwJ3CkiK4HVwN1APvB2MDJ4rdynnPP8D6zOzOesw9txROemXkcyJjapwhe3wJwXIaU5nHi/14nqhOqu1fTHGTH0QL4BDg9eHK7GabH2P2Cb33Sz3zqPAk8AzwFzcVrRnaSqeUHM4Ymych+3vL+I1Zn53HnaoTx90UDrOdqYUFk6xSk4fc+Bvy6DRm29TlQnVHek0xLwVbFcgebBCqOq1bZNVFUFJrhTTPnXtNV8MH8LV43uzp9H2kVMY0KmrASm3gFt+sGYlyA+Ii9vx6TqftObcY521hxgeX9gS1AT1UGbcgqZNGM9/52TwfDuzbntlEO9jmRM7NqxBj74CxRkwWn/tIITZtX9tj8H7heRL1R1r/8CEUkB/u6uY2pJVbnqrXms2p7HiB4tuP1UKzjGhEzuFnj9LCjeA+e9An3P9jpRnVNd0XkAOA9YLSLPAvtak/XGaWQgwIOhixf7nk9bx9Ite3hoTD8usnFxjAmtj66Cgmy4dAp0G+V1mjqpyqKjqlkiMhz4N05x2XfNRYGvgGvcjj9NLWTlFfH4tNWc1Kc15x9hoxEaEzIFO2HyJU6P0Sc/aAXHQ9WezFTVjcBpItIU6IFTeNao6q6q32kOpLCkjPQdhdwweQFxAtefcIi1UjMmVBZNhmn3Qv52pz+1o670OlGdFvAVNLfIzAlhljphe24RJz7xHXlFZQDcfXpv+rZr7HEqY2LUosnw4Xho1RdG3QKD/2QdeHrMmm2E0dz0HB78YgV7S8q5/dRDGdG9BYe1/01nD8aYYNi9Cb68FdoPhnGfQ71krxMZrOiEzcNfruSF79aRlBDHQ2P6cf7gjl5HMiZ2bV0A742DsiI462krOBHEik4YTPxuHS98t46jD2nBkxcMoHlqkteRjIlNPh/8/BxM/xuktoLLPoLWdWI0+6hhRSdESst9vPrDBqYtz2RO+i6O6NyUV8cdaQ0GjAmVkkL470WwPg0OPcMZ8TOlWbVvM+FlRSfI9hSVMn15Jq//mM6izbn0bdeI2045lHHDu1jBMSZUVn0JH14JRbth1O0w+nZrMBChrOgEyfbcIp75Zg3vzd1MSbmP9k3q89SFA/jdgEoHNDXGBMvGn+CdC6FFLzj/Veh+nNeJTBWs6ATBf2dncNdHSyn3Kecf0YELh3RiYMcmxMXZ/7SMCZnN82Dp+7DwbajfzOlloIk10Il0VnQOUlFpOX//bDmDOjXhgXP60bN1Q68jGRP79myDV04G9UHn4XDc3VZwooQVnVpSVbbmFvHn1+dSWFLO/x1/iBUcY8LBVw4L3gRfKYxPg3YDvU5kasCKTg2tzcrnxncXsnJbHiXlzlBDY4d1ZmSPFh4nMyaGlZXAsg9h7svOPTjlJdDxKGg7wOtkpoas6NTA4s27ueTFWSQmxHH5iC60bJjE4R2bcGQXa5ZpTMiowruXweqpIPEw7Bpo1Rv6nG0t1KKQFZ0AvTM7gzs+WELDpATe/stQerWxU2nGhMX8152CM+o2GH4dJNl3L5pZ0alG1p4i3p6dwZPT1zC8e3MeHtOfTs1TvI5lTOzLXA4zH4elH0CrPnD0TZBgvXlEOys6B/Dj2h1MnLGe79dk41M4qmsznr5oIC2sCxtjQi97FbxyCqBw+EVw8gNWcGKEFZ1K5BaWctkrs0lKiOOq0d0ZM6gD3Vumeh3LmNi39AP4/nHIXAJxCTD2U6dJtIkZVnT8FJaU8ebPG3l82mrKfcprlw9hSFdrJGBMyC37EGb8yyk2zbrBSQ9A7zOgaRevk5kgs6LjythZyBVvzmPFtj30bJ3KvWf0tYJjTKhlrYA5LzlTs25wzK0w/FpItoENY5UVHZwbPW+bspj0HQU8deEAzuzfzrqwMSaUMpfBl7dB+vcQn+Rctzn+PmjU1utkJsTqfNFRVe76aCk/rd/JfWf2sQ46jQklXzl8+wDMfBISkuGECTDwD9CgudfJTJjU2aJTUFzG54u3MXnuJuZt3MWoni25dGhnr2MZE5vUBxmznIKz4TsYcAkcd48d2dRBdaroKLApp5CsvCIe/GIl8zbuol3jZMYN78ItJ/eino13Y8zBU4Wd65zuanashh2rGbrue/hup3MqbeSNcMJ9Xqc0HqlTRWdTno+jH/0WcHrPePKCAfxuQDvEutIw5uCpwsrP4IenYfNsZ57EQdMu5DXsSfLpf4Sep0ByI29zGk/VqaID8Mi5/WjVKJluLRrQuXkDr+MYE93yMmHTLOeoZsWnsHMNNOnsXKvpcSI07wH1klmWlsbo/qO9TmsiQJ0qOh1S47jgyE5exzAmOqlC3jbYthi2zIU1X8O2Rc6yuARofwQcMwn6ngMJid5mNRGrThUdawVtTACK86EgCwp2wO4M2L7YKTTbl0DhDmcdiYMOQ+D4e6HrKGh9GNRL9ja3iQp1qugYYw4gP9u5QXPjD5A+E6fZjSuunjOUQK9ToE1/Z2rd167NmFqJ2qIjIlcDtwBtgWXADar6vbepjIkgqs5gZyUFUJQLGT/Bni1OgSnwm/K2Q9Fu5z1tB8CI66FlL2jQEhq2hRY97XSZCZqoLDoicgHwFHA1MNN9/FJE+qhqhqfhjAkVn88pDgU73NNf2bB1IWQtdwqL/1TqPvrKfrudpMbQoAWktoIWh0CXkU5x6XAkdBsV7p/K1DFRWXSAG4HXVPVF9/V1InIKcBVwh3exjAmQKpTudY5A9pt2V3h0p7ztvw7TXFHDtk6/ZamtILGBM9Vr8OvzfVOj9tB5hF17MZ6KuqIjIonAEcBjFRZ9DVTZB3pK4WZ48bhQRTMhNGjPHlgThdcQ1AelRVC21330m6qTUN/p+DK5MdRvCkPGQ+MOzmmvBi2cxxT3iMXuNTNRIuqKDtACiAcyK8zPBE6ouLKIjAfGuy+LZfy3S0Mbr1qNgdwI2F5N3lfdurVdXpP5LYAdVXxGuARz/1WzrT3s/8/869pkCGTdqtapzbJI3X+x+N2rbp2aLqtsXq9qPr9mVDWqJqAdTtOaYyrMvxdYVc1750ZA/kmRsL2avK+6dWu7vCbzI2HfBXv/RcK+q26d2iyL1P0Xi9+9YO+/cOy7aOxsbAdQDrSuML81sD38cWrs0wjZXk3eV926tV1e0/mRIJjZImHfVbdObZZF6v6Lxe9edevUdFnI9524lSyqiMgsYJGqjvebtxqYoqoHbEggInNVdXA4Mprgsn0X3Wz/Ra9g77tovKYD8DjwhojMBn4ArsQ57fZCNe+bFOpgJmRs30U323/RK6j7LiqPdOCXm0Nvxbk5dCnwV1Wd4W0qY4wxVYnaomOMMSb6RGNDAmOMMVHKio4xxpiwsaJTCRHpKCJpIrJcRBaLyPleZzKBE5EPRWSXiLzvdRZTNRE5Q0RWicgaEfmz13lMzdTmu2bXdCohIm2B1qq6UETaAPOAnqpa4HE0EwARGQ00BMaq6nnepjEHIiIJwHLgWJy74OcBw1V1p6fBTMBq812zI51KqOo2VV3oPt+Oc0NqM09DmYCpahqQ53UOU60hwDJV3aKq+cCXwEkeZzI1UJvvWtQVHRE5RkQ+EZEtIqIiMq6Sda4WkQ0iUiQi80Tk6IP4vCOAeFXddDC5Tfj3nQmtIOzPdsAWv9dbgPYhjm1cXn0fo67oAKk49+VcD+ytuNBvrJ0HgYHAjzhj7XTyW2ehiCytZGpXYVvNgP/wa4eh5uCEbd+ZsDjo/Wk85c3+C2ZHbuGegHxgXIV5s4AXK8xbAzxUw20nATOAy7z+OWNxCuW+c983Gnjf65+zrky12Z84Q5F86LfsSeBir3+WujgdzPexpt+1aDzSOSC/sXa+rrCo2rF2KmxHgNeAb1T1jaAFNAcUrH1nIkOA+3M2cJiItBeRVOBU4KvwpTQHEsrvY0wVHaoea6dNDbYzArgAONs9nbNQRPoFKaOpXLD2HSIyHXgPOE1ENovIsOBENDVQ7f5U1TLgJuBbYCHwL7WWa5EioO9jbb5r0drhZ0ip6kxiryDXGar6m8H8TGRS1U+AT7zOYWqnNt+1WPvDGu1j7dRltu9ii+3P6Bay/RdTRUdVS3BuMDuxwqITcVpemAhl+y622P6MbqHcf1F3es294NjDfRkHdBKRAUCOqmZQ+7F2TIjZvosttj+jm2f7z+umerVo2jca0Eqm1/zWuRpIB4pxqvUxXue2yfZdrE22P6N78mr/Wd9rxhhjwiamrukYY4yJbFZ0jDHGhI0VHWOMMWFjRccYY0zYWNExxhgTNlZ0jDHGhI0VHWOMMWFjRcdEBRF5TUQ+8+izl4rIBL/X6SJysxdZQs392dSdatS7dyXbmuC3rZj8fZmas6JjTM0dCTwfyIpeFsuD8HegLZB1kNt5zN3O5oNOZGJG1PW9ZozXVDXb6wwhlqeqB90TtKrmA/kiUh6ETCZG2JGOiUoikiQiT4pIpogUicjPIjKywjqni8gqd/kMEbnQPdXTpYrtthKRj0Vkr4hsFJE/VrLOfqfXROQKEVntfs4OEflKRBLcU3JjgdP9TjONdt/zsJttr7u9R0Uk2W+bE9zTeheKyDoRyRORj0SkRYUsY0VkiYgUu7+L1/2WNRaRSSKS5b7/OxEZXIvf9Wg3+6kiMs/N/L2IdBCRUSKySETyReQzEWle0+2busWOdEy0ehT4PfBHYD1wIzBVRA5R1W0i0gn4AHgOmAj0w+k1tzqvAZ2BE4BC4Amgy4FWdv+IP4dTXGYCTYDj3MWPAb2BZsBl7rwc97HAzb4F6IPTc28xcI/f5rvgjGB7DtAA+C/wAHCF+9lXAE8BdwKfA6n7Ptsdcv1zIBc4w/3cscA3ItJLVbcF8Luo6G/ADe423wYmA0XAeJyxV94DJgDX1WLbpq7wuqdTm2wKZMIpBp+5zxsAJcAf/JbHA+uAf7ivHwJWVNjGnTi96HY5wGf0dJeP8JvXGecP6gS/eenAze7zMTh/hBtWl7uan+9KYK3f6wk4f9Ab+827q8I6m4GHD7C944B8oH6F+QuBW6vI8cvP5jdvtPt7Odlv3rXuvEEVMi8NZJs21d3JjnRMNOoO1MMZ4wMAVS0XkZ9wjhoADgXmVHjfrGq22xvwAbP9trtRRLZW8Z5pwEZgg4h8BXwNfKCqeVV9kIich3PU0APnCCXenfxtVNVcv9dbgVbu+1sB7YH/HeAjjgBSgGznoOcXyTi/v9pY7Pc8031cUmFeq1pu29QRVnRMrAnGWB0Bb0NV80RkEHAMzqiKdwAPisiRqlppsRKRoTinyv4G/BXYDZyFczrOX2kluQK9DhuHUwSOrmTZngC3UZF/HgVQ1Yrz7DqxqZL9AzHRaB3O6bUR+2aISDwwDFjuzloJVLxoPqSa7a7E+U78sp57bahdVW9S1TJV/UZV7wD645z+O8NdXMJvj2BGAFtU9X5VnaOqa3BO4wVMVbNwrgcdf4BV5uOMZ+9T1bUVpoNtCm1MrdmRjok6qlogIv8GHhGRHcAGnCOG1vx6/8wLwI0i8hjwItAX9wI8BziSUdVVIjIVmCgi44G9OI0P9h4oi4icgXO6agbOxfpjgYbACneVdOBUEekF7MS5/rMaaC8ilwA/AScDF9Xw1wBOo4InRCQTp9FACnC8qv4LmI5z+vFjEbkVp6C2AU4Bpqvq97X4PGMOmh3pmGh1G07rqVdxLo73B05Rt1WWqm4EzsU5bbUIpyj9zX1vURXbHYdTxL4BPsVppZVexfq7gbNx/sivBG4G/uz3R/1FnAI0F8jGaaTwKfBP4Emc6yQnAvdW+xNXoKr/Bq4B/gIsBabiFFdUVYHT3J/jRWAV8C7QC+fakDGesOGqTZ0hItfj3G3fRO0ffqVEJB14VlUrXl+KqG2a6GVHOiZmicg1IjJERLqKyEU498C8ZgWnWg+4N3seVEs0EblTRPKBTkHKZWKAHemYmCUiT+DcQNoc556W/wJ/V9UST4NFMBHpjNMcHWCDqta6CxsRaYZzYyzADlXdfZDxTAywomOMMSZs7PSaMcaYsLGiY4wxJmys6BhjjAkbKzrGGGPCxoqOMcaYsLGiY4wxJmz+H0JRDWKbsOmLAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEOCAYAAACO+Hw9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAUsklEQVR4nO3de7BlZX3m8e+jXAt6EENL5zLQGhwG1BQykBECgXGGizFJmcQZY0cj4wzNBCVRSTCooTqZjBKLUYhC0jg1wemR0iqwDGJUIAo40KJNggTlEoV4g74QsaWhBdHf/LHWkdW7+pw+p8+te7/fT9Wus/d637X2+zsNz177XeuslapCktSOZyz2ACRJC8vgl6TGGPyS1BiDX5IaY/BLUmP2WOwBTMdBBx1Uy5cvX+xhSNJu4/bbb3+4qpZur223CP7ly5ezbt26xR6GJO02knx9sjaneiSpMQa/JDXG4Jekxhj8ktQYg1+SGmPwS1JjDH5JaozBL0mNMfglqTG7xV/uzsrqkyZvO+umhRuHJO0i3OOXpMYY/JLUGINfkhpj8EtSYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1BiDX5IaY/BLUmMMfklqjMEvSY0x+CWpMQa/JDXG4Jekxhj8ktQYg1+SGmPwS1JjDH5JaozBL0mNMfglqTEGvyQ1xuCXpMYY/JLUGINfkhozreBP8pNJPphkU5LvJ/lKkpMG7UmyKsmDSbYmuTHJC0a2cWCSNUk29481SZ41x/VIknZgh8Hfh/MtQICXA0cA5wAbB93OA87tlx/bt12fZMmgz5XA0cDp/eNoYM2sK5Akzcge0+hzHvBQVf32YNkDE0+SBHgTcGFVXd0vex1d+K8AVic5gi7sT6iqtX2fs4DPJTm8qu6di2IkSTs2nameVwC3JflIko1J7kjyxj7wAZ4LLAOum1ihqrYCNwPH94uOA7YAtw62ewvw2KDPNpKsTLIuybpNmzbNpCZJ0hSmE/zPA84G7gdOAy4BLgTe0Lcv639uGFlvw6BtGbCpqmqisX++cdBnG1V1eVUdU1XHLF26dBrDlCRNx3Smep4BrKuq8/vXf5/k+XTB//55G5kkaV5MZ4//IeArI8vuBg7pn6/vfx480ufgQdt6YOlgemji2MBzBn0kSQtgOsF/C3D4yLJ/BXy9f/4AXXifMtGYZB/gRJ6e018L7E831z/hOGA/tp33lyTNs+lM9bwXuDXJ24GPAC8Gfhd4G3Rz9UkuBt6W5B7gPuAddAdzr+z73J3kU3Rn+Kzst7sauNYzeiRpYe0w+Kvqi0leAbwT+CPgG/3Pywbd3g3sC1wKHAjcBpxaVY8O+qwA3gd8un99DfDGWY5fkjRD09njp6o+AXxiivYCVvWPyfo8ArxmZsOTJM01r9UjSY0x+CWpMQa/JDXG4Jekxhj8ktQYg1+SGmPwS1JjDH5JaozBL0mNMfglqTEGvyQ1xuCXpMYY/JLUGINfkhpj8EtSYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1BiDX5IaY/BLUmMMfklqjMEvSY0x+CWpMQa/JDXG4Jekxhj8ktQYg1+SGmPwS1JjDH5JaozBL0mNMfglqTEGvyQ1xuCXpMYY/JLUGINfkhpj8EtSYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1BiDX5IaY/BLUmNmHPxJzk9SSd4/WJYkq5I8mGRrkhuTvGBkvQOTrEmyuX+sSfKsOahBkjQDMwr+JC8BVgJ3jjSdB5wLnAMcC2wErk+yZNDnSuBo4PT+cTSwZueGLUnaWdMO/iQHAB8CXg88Mlge4E3AhVV1dVXdBbwOWAKs6PscQRf2K6tqbVWtBc4CfjnJ4XNUiyRpGmayx385cFVVfXZk+XOBZcB1EwuqaitwM3B8v+g4YAtw62C9W4DHBn0kSQtgj+l0SnImcBjwmu00L+t/bhhZvgH46UGfTVVVE41VVUk2DtYffc+VdNNKHHLIIdMZpiRpGna4x99PxbwTWFFVP5j/IXWq6vKqOqaqjlm6dOlCva0kjb3pTPUcBxwEfDnJU0meAk4Czu6f/3Pf7+CR9Q4G1vfP1wNL++MBwI+PDTxn0EeStACmE/wfA14EHDV4rAM+3D+/jy68T5lYIck+wIk8Pae/Ftif7kNkwnHAfmw77y9Jmmc7nOOvqu8C3x0uS/IY8J3+DB6SXAy8Lck9dB8E76A7mHtlv427k3wKWN3P3QOsBq6tqnvnpBJJ0rRM6+DuNLwb2Be4FDgQuA04taoeHfRZAbwP+HT/+hrgjXP0/pKkadqp4K+qk0deF7Cqf0y2ziNs/6wgSdIC8lo9ktQYg1+SGmPwS1JjDH5JaozBL0mNMfglqTEGvyQ1xuCXpMYY/JLUGINfkhpj8EtSYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1BiDX5IaY/BLUmMMfklqjMEvSY0x+CWpMQa/JDXG4Jekxhj8ktQYg1+SGmPwS1JjDH5JaozBL0mNMfglqTEGvyQ1xuCXpMYY/JLUGINfkhpj8EtSYwx+SWqMwS9JjTH4JakxBr8kNcbgl6TGGPyS1BiDX5IaY/BLUmMMfklqjMEvSY0x+CWpMQa/JDXG4Jekxuww+JOcn+SLSb6XZFOSjyd54UifJFmV5MEkW5PcmOQFI30OTLImyeb+sSbJs+a4HknSDkxnj/9k4DLgeOClwFPADUmePehzHnAucA5wLLARuD7JkkGfK4GjgdP7x9HAmlmOX5I0Q3vsqENVnTZ8neS1wGbgF4CPJwnwJuDCqrq67/M6uvBfAaxOcgRd2J9QVWv7PmcBn0tyeFXdO3clSZKmsjNz/Ev69R7pXz8XWAZcN9GhqrYCN9N9SwA4DtgC3DrYzi3AY4M+20iyMsm6JOs2bdq0E8OUJG3PzgT/JcAdwNr+9bL+54aRfhsGbcuATVVVE439842DPtuoqsur6piqOmbp0qU7MUxJ0vbscKpnKMl7gBPopmx+OD9DkiTNp2nv8Sd5L/Bq4KVVdf+gaX3/8+CRVQ4etK0HlvbHAya2F+A5gz6SpAUwreBPcglPh/49I80P0IX3KYP++wAn8vSc/lpgf7q5/gnHAfux7by/JGme7XCqJ8mlwGuBVwCPJJmYk99SVVuqqpJcDLwtyT3AfcA76A7mXglQVXcn+RTdGT4r+/VXA9d6Ro8kLazpzPGf3f/825Hlfwys6p+/G9gXuBQ4ELgNOLWqHh30XwG8D/h0//oa4I0zH7IkaTamcx5/ptGn6D4EVk3R5xHgNTMYmyRpHnitHklqjMEvSY0x+CWpMQa/JDXG4Jekxhj8ktQYg1+SGmPwS1JjDH5JaozBL0mNMfglqTEGvyQ1xuCXpMbM6NaLrVr+h5+YtO2fLnz5Ao5EkmbPPX5Jaox7/LM01bcB8BuBpF2Pwc+Ow1uSxolTPZLUGINfkhpj8EtSY5zjn2eeCippV+MevyQ1xuCXpMYY/JLUGINfkhpj8EtSY5o4q+fOb2/e7vJf9S92JTXIPX5JaozBL0mNMfglqTEGvyQ1xuCXpMY0cVbPrmo29wHwOj+SdpZ7/JLUGINfkhpj8EtSYwx+SWqMwS9JjTH4Jakxns65m9rRqaCe7ilpMu7xS1Jj3ONvkN8WpLa5xy9JjXGPf0zN5nIQfiOQFtDqkyZvO+umeXlL9/glqTEGvyQ1xqmeKVyz19unbP/VJ//HAo1k1zLVVJDTQNKub8GDP8nZwB8APwl8GXhTVX1uoccBOw72+Vx/qg+NVj9wvEy1tDAWNPiTvAq4BDgb+H/9z08mObKqvrGQY9mdzfYDa2dN5wNnNuG9q/Jgt8bNQu/xvwW4oqo+0L8+J8npwO8A5y/wWBbVYoX3bMx2zIv2TWWqsyYAzrppcT6wpjEuaT4sWPAn2Qv4N8BFI03XAccv1Di0eObzw275H071vpunXvmCo7hmr51/7zsv2Pm6fu6nD5i8cUcfDDtw57cnr/vn/uSOnd/uBUdN2T5lTfNpPj8o5/GUyyn/nWa15cmlquZp0yNvlPwU8G3gpKq6ebD8AuC3qurwkf4rgZX9y8OBe3fyrQ8CHt7JdXdX1jz+WqsXrHmmDq2qpdtr2GXP6qmqy4HLZ7udJOuq6pg5GNJuw5rHX2v1gjXPpYU8j/9h4IfAwSPLDwbWL+A4JKlpCxb8VfUkcDtwykjTKcCtCzUOSWrdQk/1vAdYk+QLwC3AfwN+CvjLeXzPWU8X7Yasefy1Vi9Y85xZsIO7P37D7g+4zqP7A667gDcPD/ZKkubXgge/JGlxeZE2SWqMwS9JjRnb4E9ydpIHknw/ye1JTlzsMc2VJL+Y5Jok305SSc4YaU+SVUkeTLI1yY1JXrBIw50TSc5P8sUk30uyKcnHk7xwpM9Y1Z3kDUnu7Gv+XpK1SV4+aB+rekf1/+aV5P2DZWNVc19LjTzWD9rnpd6xDP7BxeDeCbyY7nTRTyY5ZFEHNnf2pzsw/nvA1u20nwecC5wDHAtsBK5PsmTBRjj3TgYuo7u8x0uBp4Abkjx70Gfc6v4W8FbgaOAY4DPAx5JM/CX/uNX7Y0leQveX+3eONI1jzffSnewy8XjRoG1+6q2qsXsAtwEfGFn2j8C7Fnts81DrFuCMwesADwFvHyzbF3gUOGuxxzuHde9P9weBv9JY3d8BzhrneoEDgK8B/w64EXj/uP4bA6uAuyZpm7d6x26Pf3AxuOtGmlq5GNxzgWUM6q+qrcDNjFf9S+i+sT7Svx7rupM8M8lv0n3g3cp413s5cFVVfXZk+bjW/Lx+KueBJB9O8rx++bzVO3bBT3dRo2cCG0aWb6D7JY67iRrHvf5LgDuAtf3rsaw7yYuSbAGeoPtDx1+rqn9gfOs9EzgMeMd2msex5tuAM4DTgTPp6rg1yU8wj/XushdpkyaT5D3ACcAJVfXDxR7PPLsXOIpu+uOVwAeTnLyI45k3SQ6nOy53QlX9YLHHsxCq6pPD10k+D9wPvA74/Hy97zju8bd+MbiJGsey/iTvBV4NvLSq7h80jWXdVfVkVX21qm6vqvPpvuW8mfGs9zi6b+xfTvJUkqeAk4Cz++f/3Pcbp5q3UVVb6G5J+3zm8d947IK/vBjcA3T/Ufy4/iT7ACeym9ef5BKeDv17RprHtu4RzwD2Zjzr/RjdGS1HDR7rgA/3z+9j/GreRl/Pv6Y7qDt//8aLfVR7no6Uvwp4EvivwBF088Fb6G5MsOjjm4P69ufp/zEeBy7onx/St78V2Az8OvBCuv9xHgSWLPbYZ1HzpcD36E7lXDZ47D/oM1Z1Axf2/5MvpwvEdwE/Al42jvVO8ju4kf6snnGsme6OhCfRHcj9t8C1/X/nh85nvYte+Dz+Qs8G/onuoNjtwC8u9pjmsLaTgdrO44q+PXSniT0EfB+4CXjhYo97ljVvr94CVg36jFXdwBXA1/v/hjcCNwCnjWu9k/wORoN/rGoeBPmTdHcovBo4cr7r9SJtktSYsZvjlyRNzeCXpMYY/JLUGINfkhpj8EtSYwx+SWqMwS/thgY37fj+HGzrisH2XjkX49OuzeDXghkJmB8kuT/JRUn2W+yx7abOBA6dg+38Ht0NQNQIr86phXYD8FpgT7rLEfwvYD/gdxZzUDORZM8auXpkkr2qu07UTLe1U+v1vltVo5fsnbGq2gxsTjLbTWk34R6/FtoTVbW+qr5ZVVcCHwJeAZBk7yQXJ9mQ7l7Jn09ywsSKSfZM8uf9TSueSPLNJBcO2vdK8mdJvpXk8f4evadNNZj+nqbnJflaf0/Tf0jymkH78v4byquTfCbJVuCs/tvLtUnemuRbdLdJnLh+/g39tr7T9ztgsL3J1vv1dPfXnVjvpiSjV2WcUpIzkmxJ8rIk9/S/g2uSHJDklUn+McnmJGuS7DuTbWu8uMevxbaVbu8f4N3AfwJeT3dN8rcAn0ry/Kp6CPhd4NeA36S7DtPPAIcPtvVXwM8CK+gC9ZeAjyc5tqq+NMn7/yndde7fQHft++OADyR5pKo+Mej3LuD3gf8C/IDuHrgn0V1A63S6z5D9gE8DXwB+Hng28AHgfwO/MdjW6HrL6K7Zcj7dtVr2B16yg9/bZPamu0frbwF79du7mu73/BvATwAfpbuW1f/cyffQ7m6xL1Lko50H3UXHrh28/nm6+yd8hG6650ngtwftz6S79+qf9q//HPhb6K4xNbLtn6W7cuUhI8s/Blw2yXj2owvEE0eWXwz8Tf98Od3F4M7dTi2bgL0Hy86kC/Qlg2Un9+sfNsV6R/d9Dp3B77KAV44sO6Nffvhg2UV096c4aLJ/h6m26WM8H+7xa6Gd3t9KcA+6Pf2/Bs6hC+49gVsmOlbVD5OsBY7sF10BXA/cl+Q64G+AT1bVj+jCM8BXRuaq9wY+M8lYjgT2oftWMbxa4Z503yiG1m1n/buq6onB6yOAO6vq0cGyW+k+kI4EvjrJel+iO/ZxV1/XDXT3nN00ybin8kRV3Tt4vQFYX1UPjyw7EjXL4NdCuxlYSTdd8mD1B0n76Y7JdLujVX+XZDlwGvDvgQ8CX0pyCt3xqgKO7bc9tHWS7U4c4/oV4BsjbaPbeGw7629v2WSGHyzbrNd/wJ1KN71zKt100ruSnFSTT1FN5qntvO9oLYXH95pm8GuhPV5VX93O8q/RTfX8Qv+cJM+km3O/cqJTvzd9FXBVkivo7kt6GPD3dHv8y6rqs9Mcy1fornV/aFVN9q1gJu4GXp9kyWCv/3i6kL17qhWrquhuHL82yZ/Q3X7vVXTfBqQ5ZfBrl1BVjyX5C+DPkjxMd9u5N9PdX/QygCRvobshxR10e7Er6O5W9K2qejzJh4ArkpwL/B3dwdWTgfur6qPbec9Hk1wEXJRufuhmnj6w+qOqunyGZXwI+GPg/yS5ADgQWA18dJIPO/q6XgL8B7oDwxuAFwP/ku6DSZpzBr92JW/tf/4V8Cy6vfjTqzujB+BR4A/obkRdffvLqurxvv0/A2+nOzvoZ4Dv0J1hM9U3gD+iC9vfB/6C7oPkjn4bM9J/+JxGd3D4C3R3TPpruj+Qmspmum8659DV/U3gv1fV/53pGKTp8A5c0m6oPxj9H6vqql15m9o1eYBH2n2t6afFZiXJX/ZnWqkR7vFLu6Ekh/VPf1RV989yW88B/kX/8qGqmsnZStoNGfyS1BineiSpMQa/JDXG4Jekxhj8ktQYg1+SGvP/AUT8jvVJnUWjAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "x, y = cumulative_recall(err_pred)\n",
    "plt.plot(x, y)\n",
    "recall = y[np.searchsorted(x, [0.25, 0.5, 5])]\n",
    "print('Prediction:', recall)\n",
    "\n",
    "x, y = cumulative_recall(err_init)\n",
    "plt.plot(x, y)\n",
    "recall = y[np.searchsorted(x, [0.25, 0.5, 5])]\n",
    "print('Initial:   ', recall)\n",
    "\n",
    "plt.xlim([0.01, 10])\n",
    "plt.ylim([0, 100])\n",
    "plt.grid()\n",
    "plt.xscale('log')\n",
    "plt.xlabel('log distance [m]')\n",
    "plt.ylabel('Cumulative recall [%]');\n",
    "plt.title(f'Aachen');\n",
    "\n",
    "plt.figure()\n",
    "plt.hist(np.clip(err_init, a_max=50, a_min=None), bins=40, alpha=1.);\n",
    "plt.hist(np.clip(err_pred, a_max=50, a_min=None), bins=40, alpha=.8);\n",
    "plt.xlabel('Pose errors [m]');"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
